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Abstract 

The venerable theory of electrokinetic phenomena rests on the hypothesis of a dilute solution 
of point-like ions in quasi-equilibrium with a weakly charged surface, whose potential relative 
to the bulk is of order the thermal voltage {kT le ^ 25 mV at room temperature). In nonlinear 
electrokinetic phenomena, such as AC or induced-charge electro-osmosis (ACEO, ICEO) and 
induced-charge electrophoresis (ICEP), several Volts ^ lOOkT je are applied to polarizable sur- 
faces in microscopic geometries, and the resulting electric fields and induced surface charges are 
large enough to violate the assumptions of the classical theory. In this article, we review the 
experimental and theoretical literatures, highlight discrepancies between theory and experiment, 
introduce possible modifications of the theory, and analyze their consequences. We argue that, 
in response to a large applied voltage, the "compact layer" and "shear plane" eff'ectively advance 
into the liquid, due to the crowding of counter-ions. Using simple continuum models, we pre- 
dict two general trends at large voltages: (i) ionic crowding against a blocking surface expands 
the diff'use double layer and thus decreases its diff'erential capacitance, and (ii) a charge-induced 
viscosity increase near the surface reduces the electro-osmotic mobility; each trend is enhanced 
by dielectric saturation. The first effect is able to predict high-frequency flow reversal in ACEO 
pumps, while the second may explain the decay of ICEO flow with increasing salt concentration. 
Through several colloidal examples, such as ICEP of an uncharged metal sphere in an asym- 
metric electrolyte, we show that nonlinear electrokinetic phenomena are generally ion-specific. 
Similar theoretical issues arise in nanofluidics (due to confinement) and ionic liquids (due to 
the lack of solvent), so the paper concludes with a general framework of modified electrokinetic 
equations for finite- sized ions. 
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1. Introduction 



1.1. Nonlinear ''induced-charge" electrokinetic phenomena 

Due to favorable scaling with miniaturization, electrokinetic phenomena are finding many new 
applications in microfluidics CJEIEI. but often in new situations that raise fundamental theo- 
retical questions. The classical theory of electrokinetics, dating back to Helmholtz and Smolu- 
chowski a century ago |4|, was developed for the eff'ective linear hydrodynamic slip driven by an 
electric field past a surface in chemical equilibrium with the solution, whose double-layer voltage 
is of order the thermal voltage, kT je = 25 mV, and approximately constant [ 5^/61171 [8l[9l[TQll. 
The discovery of AC electro-osmotic flow (ACEO) over micro-electrodes 1 1 1 , 12 , 13l has shifted 
attention to a new nonlinear regime lfT4l . where the induced double-layer voltage is typically sev- 
eral Volts ^ 100 kT/e, oscillating at frequencies up to 100 kHz, and nonuniform at the micron 
scale. Related phenomena of induced-charge electro-osmosis (ICEO) 1X51 HSl [EJ also occur 
around polarizable particles |[l8l[T9l and microstructures |[20l|2Tl (in AC or DC fields), as well 
as driven biological membranes f22l . Due to broken symmetries in ICEO flow, asymmetric col- 
loidal particles undergo nonlinear, induced-charge electrophoresis (ICEP) 1 15 , 23 , 24, 25 1. Some 
of these fundamental nonlinear electrokinetic phenomena are illustrated in Fig. [T] 

A "Standard Model" (outlined below) has emerged to describe a wide variety of induced- 
charge electrokinetic phenomena, but some crucial aspects remain unexplained. In their pioneer- 
ing work 25 years ago in the USSR, which went unnoticed in the West until recently [|T5l[T6l , 
V. A. Murtsovkin, A. S. Dukhin and collaborators first predicted quadrupolar flow (which we 
call "ICEO") around a polarizable sphere in a uniform electric field |28| and observed the 
phenomenon using mercury drops 1291 and metal particles |30|, although the flow was some- 
times in the opposite direction to the theory, as discussed below. (See Ref. fTSl for a re- 
view.) More recently, in microfluidics, Ramos et al. observed and modeled ACEO flows over 
a pair of micro-electrodes, and the theory over-predicted the observed velocity by an order of 
magnitude (TT] [31] [32j [33l. Around the same time, Ajdari used a similar model to predict 
ACEO pumping by asymmetric electrode arrays 1 12L which was demonstrated using planar 
electrodes of unequal widths and gaps O |35] [36| [371 [38j [39l, but the model cannot predict 
experimentally observed flow reversal at high frequency and loss of flow at high salt concentra- 
tion 1 38, 40, 41], even if extended to large voltages |42l|43l. The same model has also been 
used to predict faster three-dimensional ACEO pump geometries 1441 . prior to experimental 
verification ||40l |45] |46l |47l |48l, but again the data departs from the theory at large voltages. 
Discrepancies between theory and experiments, including flow reversal, also arise in traveling- 
wave electro-osmosis (TWEO) for electrode arrays applying a wave-like four-phase voltage 
pulse ||49j|50l[5TJ[52l. Recent observations of ICEO flow around metal microstructures 12011211 . 
ICEP rotation of metal rods | 53 1, ICEP translation of metallo-dielectric particles | 25 1 have like- 
wise confirmed qualitative theoretical predictions |[T5l[T6l|23]|54j|55l, while exhibiting the same 
poorly understood decay of the velocity at high salt concentration. We conclude that there are 
fundamental gaps in our understanding of induced-charge electrokinetic phenomena. 

In this article, we review recent experimental and theoretical work in this growing area of 
nonlinear electrokinetics, as well as some possibly relevant literatures from other fields, and we 
consider a number of possible modifications of classical electrokinetic theory. Some of these 
ideas are new, while others were proposed long ago by O. Stern, J. J. Bikerman, J. Lyklema, 
and others, and eff'ectively forgotten. We build the case that at least some failures of the Stan- 
dard Model can be attributed to the breakdown of the dilute-solution approximation at large 
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Figure 1: Examples of nonlinear electrokinetic phenomena, driven by induced charge (+, -)in the diffuse part of the 
electrochemical double layer at polarizable, blocking surfaces, subject to an applied electric field E or voltage V. (a) 
Induced-charge electro-osmosis (ICEO) around a metal post 1 15][l6][20][E k (b) induced-charge electrophoresis (ICEP) 
of a metal/insulator Janus particle | 23 25 1, (c) a nonlinear electrokinetic jet of ICEO flow at a sharp corner in a dielectric 
microchannel t26|i27J, and (d) AC electro-osmosis (ACEO) over a pair of microelectrodes 1 1 1 121. 
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induced voltages. Using simple models, we predict two general effects associated with coun- 
terion crowding - decay of the double-layer capacitance and reduction of the electro-osmotic 
mobility - which begin to explain the experimental data. Our models, although incomplete, also 
imply generic new ion- specific nonlinear electrokinetic phenomena at large voltages related to 
atomic-level details of polarizable solid/electrolyte interfaces. 



7.2. Scope and context of the article 

We first presented these ideas in a paper at the ELKIN International Electrokinetics Sympo- 
sium in Nancy, France in June 2006 [56] and in a letter, which was archived online in March 
2007 271 and recently published |[58l . The present article is a review article with original 
material, built around the letter, where its basic arguments are further developed as follows: 

1 . We begin with a critical review of recent studies of induce-charge electrokinetic phenomena 
in section |2] By compiling results from the literature and performing our own simulations 
of other experiments, we systematically compare theory and experiment across a wide range 
of nonlinear electrokinetic phenomena. To motivate modified electrokinetic models, we also 
review various concentrated- solution theories from electrochemistry and electrokinetics in 
sections 3 and 4. 

2. In our original letter, the theoretical predictions of steric eff'ects of finite ion sizes in elec- 
trokinetics were based on what we call the "Bikerman's model" below [59l[60l, a simple 
lattice-gas approach that allows analytical results. Here, we develop a general, mean-field 
local-density theory of volume constraints and illustrate it with hard-sphere liquid mod- 
els [61, ,62 1 . In addition to the charge-induced thickening eff'ect from the original letter, 
we also consider the field-induced viscoelectric eff'ect in the solvent proposed by Lyklema 
and Overbeek (631 IMI. in conjunction with our models for steric eff'ects. We also consider 
dielectric relaxation of the solution, which tends to enhance these eff'ects. 

3. We provide additional examples of new electrokinetic phenomena predicted by our models 
at large voltages. In the letter |57|, we predicted high-frequency flow reversal in ACEO 
(Fig. [To] below) and decay of ICEO flow at high concentration (Fig. [17]). Here, we also 
predict two mechanisms for ion- specific, field-dependent mobility of polarizable colloids 
at large voltages. The first has to do with crowding eff'ects on the redistribution of double- 
layer charge due to nonlinear capacitance, as noted by A. S. Dukhin f65",^6 | (Fig. [12]). The 



second results from a novel ion-specific viscosity increase at high charge density (Fig.[2Q|. 
4. We present a general theoretical framework of modified electrokinetic equations and bound- 
ary conditions for concentrated solutions and/or large voltages in section[5] which could find 
many other applications in nonlinear electrochemical relaxation or electrokinetics. . 

In spite of these major changes, the goal of the paper remains the same: to provide an overview of 
various physical aspects of electrokinetic phenomena, not captured by classical theories, which 
become important at large induced voltages. Here, we focus on general concepts, mathematical 
models, and simple analytical predictions. Detailed studies of some particular phenomena will 
appear elsewhere, e.g. Ref. 1671 on high-frequency flow reversal in AC electro-osmosis. 

There have been a few other attempts to go beyond dilute solution theory in electrokinetics, 
but in the rather difl'erent context of linear "fixed-charge" flows in nanochannels at low surface 
potentials. The first electrokinetic theories of this type may be those of Cervera et al. I68ll69l , 
who used Bikerman's modified Poisson-Boltzmann (MPB) theory to account for the crowding 
of finite-sized ions during transport by conduction and electro-osmosis through a membrane 

4 



nanopore. Independently, J. J. Horno et al. ||70l |7T] |72j |73l also used Bikerman's model (albeit, 
attributed to others |[74j [75l |76l - see below) to analyze linear electrophoresis of colloids in a 
concentrated electrolyte. Recently, Liu et al. fTTl numerically implemented a more complicated 
MPB theory fTSl |79l [80l to predict effects of finite ion sizes, electrostatic correlations, and di- 
electric image forces on electro-osmotic flow and streaming potential in a nanochannel. In these 
studies of linear electrokinetic phenomena, efl'ects beyond the dilute solution approximation can 
arise due to nano-confinement, but, as we first noted in Ref. | 60'1, much stronger and possibly dif- 
ferent crowding eff'ects also arise due to large induced voltages, regardless of confinement or bulk 
salt concentration. Our goal here is to make a crude first attempt to understand the implications of 
ion crowding for nonlinear electrokinetic phenomena, using simple mean-field approximations 
that permit analytical solutions in the thin double-layer limit. 

Similar models for double-layer charging dynamics are also starting to be developed for ionic 
liquids and molten salts I8l]|82j|83]|84|, since describing ion crowding is paramount in the ab- 
sence of a solvent. Kornyshev recently suggested using what we call the "Bikerman-Freise" (BF) 
mean-field theory below to describe the diff'erential capacitance of the double layer, with the bulk 
volume fraction of ions appearing as a fitting parameter to allow for a slightly diff'erent density 
of ions | 81 1. (An equivalent lattice-gas model was also developed long ago for the double layer 
in an ionic crystal by Grimley and Mott |[85l[86||87l, and a complete history of related models is 
given below in Section [3^1.2[ ) The BF capacitance formula, extended to allow for a thin dielectric 
Stern layer, has managed to fit recent experiments and simulations of simple ionic liquids rather 
well, especially at large voltages ||82j|83l[88l. However, we are not aware of any work addressing 
electrokinetic phenomena in ionic liquids, so perhaps the mean-field electro-hydrodynamic mod- 
els developed here for concentrated electrolytes at large voltages might provide a useful starting 
point, in the limit of nearly close packing of ions. 

As a by-product of this work, our attempts to model nanoscale phenomena in nonlinear elec- 
trokinetics may also have broader applicability in nanofluidics | 89 , 90J. Dilute- solution theory 
remains the state-of-the-art in mathematical modeling, and the main focus of the field so far 
has been on eff'ects of geometrical confinement, especially with overlapping double layers. The 
classical Poisson-Nemst-Planck and Navier-Stokes equations provide a useful first approxima- 
tion to understand nanochannel transport phenomena, such as charge selectivity 19111921 19311941 , 
mechanical-to-electrical power conversion efficiency 1*951 [96l [97l [98l l99l IIOOL current- voltage 
characteristics |101|, and nonlinear ion-profile dynamics 11021 11Q3L but in some cases it may 
be essential to introduce new physics into the governing equations and boundary conditions to 
account for crowding effects and strong surface interactions. Molecular dynamics simulations 
of nanochannel electrokinetics provide crucial insights and can be used to test and guide the 
development of modified continuum models |[T0l|T05l|T06l[l07l[l0^ . 

We stress that there are other important, developing areas of nonlinear electrokinetics, which 
are related, but outside the scope of this article. For example, we do not discuss the non- 
linear electrophoresis of fixed-charge particles IllQI lllll 1661 , which can result from surface 
conduction fTTT, TTSl [mj [TTSl (Du > 0) and advection-diff'usion fTTOl (Pe > 0), although 
we will analyze a diff'erent mechanism for field-dependent electrophoretic mobility for po- 
larizable particles 1651 . We also neglect nonlinear electrokinetic phenomena associated with 
strong salt concentration gradients, such as second-kind "superfast" electrophoresis 11161 11171 
fTTSl ITT9I fT20l [ml . electro-osmotic ffuid instability 1122 [Ell \lM [123 [Ml, and 
concentration-polarization shocks 11021 11Q3L which can now be directly observed in microflu- 
idic systems ri25'J 126lllQllllQ3l and porous bead packings 1127111281 . These efl'ects result from 
"super-limiting" normal current into a polarized surface, membrane, or nanochannel, which de- 
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pletes the local salt concentration significantly and forces the diff'use charge out of equilibrium. 
In such cases, dilute solution theory is likely to remain valid up to large applied voltages in the 
regions of low ionic strength away from the surface, where the flow is mostly generated. In 
contrast, our focus here is on metal structures and electrodes that do not sustain normal current, 
e.g. due to AC forcing and insufficient voltage to trigger Faradaic reactions. As a result, the 
applied voltage leads mostly to capacitive charging of the diff'use layer, and thus potentially to 
the crowding of counter-ions attracted to the surface. 

The article is organized as follows. In section [2j we review the standard low- voltage model 
for induced-charge electrokinetic phenomena and its failure to explain certain key experimental 
trends. We then review various attempts to go beyond dilute solution theory in electrochemistry 
and electrokinetics and analyze the eff'ects of two types of new physics in nonlinear electroki- 
netic phenomena at large voltages: In section [3j we build on our recent work on diff'use-charge 
dynamics at large applied voltages ll29|[6Ql [T3QI to argue that the crowding of counterions plays 
a major role in induced-charge electrokinetic phenomena by reducing the double-layer capaci- 
tance in ways that are ion-specific and concentration-dependent; In section |4j we postulate that 
the local viscosity of the solution grows with increasing charge density, which in turn decreases 
the electro-osmotic mobility at high voltage and/or concentration and introduces another source 
of ion specificity. Finally, in section |5j we present a theoretical framework of modified elec- 
trokinetic equations, which underlies the results in sections [3]|4] and can be applied to general 
situations. 

2. Background: Theory versus Experiment 

2. 1 . The Standard Model 

We begin by summarizing the "Standard Model" of induced-charge / AC electrokinetics, used 
in most theoretical studies, and bringing out some crucial experimental trends it fails to capture. 
The general starting point for the Standard Model is the coupling of the Poisson-Nemst-Planck 
(PNP) equations of ion transport to the Navier Stokes equations of viscous fluid flow. ICEO 
flows are rather complex, so many simplifications from this starting point have been made to 
arrive at an operational model 1 18^, jT| [T2| O [23J. For thin double layers (DL) compared to 
the geometrical length scales, the Standard Model is based on the assumption of "linear" or 
"weakly nonlinear" charging dynamics [1291, which further requires that the applied voltage is 
small enough not to significantly perturb the bulk salt concentration, whether by double-layer 
salt adsorption or Faradaic reaction currents. In this regime, the problem is greatly simplified, 
and the electrokinetic problem decouples into one of electrochemical relaxation and another of 
viscous flow: 

1. Electrochemical relaxation. - The first step is to solve Laplace's equation for the electro- 
static potential across the bulk resistance, 

V • J = V • (o-E) = -crV^^ = (1) 

assuming Ohm's Law with a constant conductivity cr. A capacitance-like boundary con- 
dition for charging of the double layer at a blocking surface (which cannot pass normal 
current) then closes the "RC circuit" I129L 



(2) 



where the local diffuse-layer voltage drop (surface minus bulk) responds to the nor- 

mal electric field En = -h- V0. In the Standard Model, the bulk conductivity cr and diff'use- 
layer capacitance Cd are usually taken to be constants, although these assumptions can be 
relaxed f42l |43l 11311 . The diff'use layer capacitance is calculated from the PNP equations 
by applying the justifiable assumption that the thin double layers are in thermal equilibrium; 
see section 3.1. A compact Stem layer or dielectric surface coating of constant capacitance 
Cs is often included 112, 129ll42l. so that only part of the total double-layer voltage A0 is 
dropped across the diff'use layer "capacitor", 

where S = Cd/Cs is the diff'use-layer to compact-layer capacitance ratio. 
2. Viscous flow. - The second step is to solve for a (possibly unsteady) Stokes flow, 

P'n^= -^P + V • U = 0, (4) 

ot 

with the Helmholtz-Smoluchowski (HS) boundary condition for efl'ective fluid slip outside 
the double layer, 

u, = -bE, = -^-^E, (5) 

where is the tangential field, b = Sb^/rjiy is the electro-osmotic mobility, f is the zeta 
potential at the shear plane (= in the simplest models), and Sh, r]iy, and pm are the per- 
mittivity, viscosity, and mass density of the bulk solvent. Osmotic pressure gradients, which 
would modify the slip formula 1123111241 . are neglected since the bulk salt concentration is 
assumed to be uniform. 

Although this model can be rigorously justified only for very small voltages, ^ kT/e, in a. 
dilute solution |[32l[T6l[T29]| , it manages to describe many features of ICEO flows at much larger 
voltages. 

There has been extensive theoretical work using the Standard Model, and it provides the ba- 
sis for most of our understanding of induced-charge electrokinetics. In recent years, it has been 
widely used to model nonlinear electrokinetic phenomena in microfluidic devices, such as ACEO 
flows around electrode pairs CU HI |32l |33 II3T| and arrays (121 |3l |44| |45| |46| gTl \lM 11341 . 
traveling- wave electro-osmotic flows (TWEO) 149L .5()_, 52J, ICEO interactions between dielectric 
particles and electrodes |[l35l[l36l[l37l[l38l, ICEO flow around metal structures |[T5l[T6ll20lfT39l 
[1401 [Hi] [111 and dielectric corners f26','2n and particles |[l6l[l44l[l45l, fixed-potential 
ICEO around electrodes with a DC bias 1 16, 146], ICEP motion of polarizable asymmetric par- 
ticles (HI [mm [231, collections of interacting particles |[53l[54l[Tt7l[T48l[T49l . particles near 
walls 11501 [55l. and particles in field gradients |23 1. The Standard Model has had many suc- 
cesses in describing all of these phenomena, but it also has some fundamental shortcomings, 
when compared to experimental data. 

One possibility is that the underlying PNP/NS equations and boundary conditions are physi- 
cally accurate, but the thin-DL approximation introduces large mathematical errors. A number 
of theoretical studies have allowed for arbitrary DL thickness in a dilute solution while solving 
the linearized equations of ion transport and fluid flow in the regime of low voltages. This mod- 
eling approach has been applied to ACEO [32] and TWEO 115 1 J flows over electrode arrays 
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and ICEP particle motion in uniform 11521 or non-uniform flSSI fields. For the model problem 
of ICEO flow around a thick metal stripe on a flat wall in a parallel electric field 11291 11461 , 
Gregersen et al. II143II have recently compared full numerical solutions of the Poisson-Nemst- 
Planck (PNP) and Navier- Stokes (NS) equations with thin double-layer approximations in both 
the linear regime (Standard Model) and the "weakly nonlinear" regime 1 129, .1311 (including 
tangential surface conduction, but not diff'use-layer salt adsorption, surface reactions, and bulk 
concentration polarization, discussed below). For their model problem, they concluded that for 
micron- scale electrodes, the (outer) boundary-layer approximations can over-estimate ICEO ve- 
locities by 10% for thin double layers, and by 100% for double-layer thickness comparable to the 
electrode height. These errors could be reduced by constructing uniformly valid approximations 
for finite double layer thickness 1 154.. 44 1 (adding double layer contributions and subtracting the 
overlap), but the main point for us is that mathematical errors in thin double layer approxima- 
tions (compared to full PNP/NS numerical solutions) vanish in the thin double layer limit and 
are small (of order Ad/L ^ 1) for typical experimental situations in microfluidics. In particular, 
we cannot attribute the systematic and large (sometimes order of magnitude) discrepancies be- 
tween theory and experiment discussed below to the thin double-layer approximation, especially 
if the asymptotic analysis is done carefully, going beyond the leading-order low-voltage approx- 
imation of the Standard Model. Instead, in this paper, we build the case that at least some of 
the discrepancies may be attributable to the breakdown of the underlying PBP/NS equations of 
dilute solution theory (and its boundary conditions), close to a highly charged surface. 

2.2. Open questions 

2.2.1. The correction factor'' 

Low-voltage, dilute- solution theories in nonlinear electrokinetics tend to over-predict fluid 
velocities, compared to experiments. A crude way to quantify this efl'ect in the Standard Model 
is to multiply the HS slip velocity ([5} on all surfaces by a fitting parameter A, the "correction 
factor" introduced by Green et al. |3T, "331. This approach works best at low voltages and in 
very dilute solutions, but even in such a regime, we should stress that it is generally impossible to 
choose A to fit complete flow profiles or multiple experimental trends, e.g. velocity vs. voltage 
and frequency, at the same time. Nevertheless, one can often make a meaningful fit of A to 
reproduce a key quantity, such as the maximum flow rate or particle velocity. Such a quantitative 
test of the model has been attempted for a number of data sets EH [33l [20l [l32l [1461 15711251 l55l . 
but there has been no attempt to synthesize results from difl'erent types of experiments to seek 
general trends in the correction factor. 

As a background for our study, we provide a critical evaluation of the Standard Model based on 
A values for a wide range of experimental situations. In Table [T] we have complied all available 
results from the literature. We have also added many entries by fitting our own Standard-Model 
simulations to published data, for which no comparison has previously been done. 

It is striking that A is never larger than unity and can be orders of magnitude smaller. We 
managed to find only one published measurement where the Standard Model correctly predicts 
the maximum of the observed flow (A = 1), from a recent experiment on ACEO pumping of 
micromolar KCl by a planar, gold electrode array at relatively low voltage | 41 1, but even in that 
data set the model fails to predict weak flow reversal at high frequency and salt concentration 
dependence (see below). Remarkably, there has not yet been a single ICEO experiment where 
the model has been able to predict, or even to fit, how the velocity depends on the basic operating 
conditions - voltage, AC frequency, and salt concentration - let alone the dependence on surface 
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Table 1 : Comparison of the standard low-voltage model of induced-charge electrokinetic phenomena with experimental 
data (column 1) for a wide range of situations (column 2), although limited to a small set of aqueous electrolytes (column 
3) at low bulk salt concentrations cq (column 4). In each case, the nominal maximum induced double-layer voltage 
Vmax is estimated (column 5). A crude comparison with the Standard Model is made by multiplying the predicted slip 
velocity j5} everywhere by a constant factor A (column 6) for a given cq and Ymax- In addition to A values from the cited 
papers, we have added entries to the table, indicated by *, by fitting our own standard-model simulations to published 
experimental data. Estimates indicated by * assume a frequency-dependent constant-phase-angle impedance for the 
double layer, and those labeled by ^ are affected by particle-wall interactions, which are not fully understood. In each 
case, we also estimate the maximum induced zeta potential ^^ajc = ^max^^ in Volts (column 7) and in units of thermal 
voltage kTfe (column 8). 
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and bulk chemistries. The greatest discrepancies come from ACEO pumping by a disk-annulus 
electrode pair 1 132] (A = 0.0025) and fixed-potential ICEO around a metal stripe [1461 (A = 
0.005), both in millimolar KCl and at high induced voltages Vmax- 

In Table [T] we have also used A to convert the maximum nominal voltage Vmax induced 
across the double layer in each experiment to a maximum zeta potential ^max = ^Vmax- The 
range of ^max is much smaller for A, but still quite significant. It is clear that ^max rarely exceeds 
lOkT/e, regardless of the applied voltage. For very dilute solutions, the largest value in the 
Table, ^max = 0.75V = 30kT/e, comes from ACEO pumping of micromolar KCl |41 1, while 
the smallest values, ^^ax < 0.5kT/e, come from ICEP of gold-latex Janus particles in millimolar 
NaCl. 

The values of A and ^max f^om all the diff'erent experimental situations in Table [T] are plotted 
versus cq andV;^^;^ in Fig. [2] and some general trends become evident. In Fig. |2jb), we see that 
^max decays strongly with increasing salt concentration and becomes negligible in most exper- 
iments above 10 mM. In Fig. [2jc), we see that ^max exhibits sub-linear growth with Vmax and 
appears to saturate below ten times times the thermal voltage (lOkT/e = 0.25 V at room temper- 
ature), or much lower values at high salt concentration. Even with applied voltages up to 10 Volts 
in dilute solutions, the eff'ective maximum zeta potential tends to stay well below 1 Volt. From 
the perspective of classical electrokinetic theory, this implies that most of the voltage applied to 
the double layer is dropped across the immobile, inner "compact layer", rather than the mobile 
outer "diff'use layer", where electro-osmotic flow is generated. 

This efl'ect can be qualitatively, but not quantitatively, understood using the Standard Model. 
Many authors have assumed a uniform, uncharged Stern layer (or dielectric thin film) of permit- 
tivity Es and thickness hs = ss/Cs, acting as a capacitor in series with the diff'use layer. Via 
Eq. ([3]), this model implies 

A = , with S = — = = — , (6) 

1 -\-S Cs ss Ad Ad 

where is the Debye-Htickel screening length (difl'use-layer thickness), which takes the form 



^ 2{zeycQ 

for a z : z electrolyte, and As is an efl'ective width for the Stern layer, if it were a capacitor 
with the same dielectric constant as the bulk. Inclusion of the Stern layer only transfers the 
large, unexplained variation in the correction factor A to the parameter As (or Cs = Sb/As) 
without any theoretical prediction of why it should vary so much with voltage, concentration, 
and geometry. Using these kinds of equivalent circuit models applied to difl'erential capacitance 
measurements 1 155], electrochemists sometimes infer a tenfold reduction in permittivity in the 
Stern layer versus bulk water, st/ss ~ 10, but, even if this were always true, it would still be 
hard to explain the data. For many experimental situations in Table [T] the screening length Ad 
is tens of nanometers, or hundreds of molecular widths, and the efl'ective Stern-layer width As 
would need to be much larger - up to several microns - to predict the observed values of A ^ 1 . 
In contrast, if we take the physical picture of a Stern monolayer literally, then hs should be only 
a few Angstroms, and As at most a few nanometers, so there is no way to justify the model. As 
noted in early papers by Brown et al. 1341 and Green et al. |31 1 , it is clear that the efl'ective 
difl'use-layer voltage (or induced zeta potential) is not properly described by the Standard Model 
under typical experimental conditions. 
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Figure 2: General trends in the under-prediction of ICEO flow velocity by the Standard Model from Table ^ compared 
with (solid and dashed) scaling curves, simply to guide the eye. (a) Log-log plot of A vs. cq compared with the curves 
A = ^JlO~^mM/co (dashed) and = lnio(10mM/co)/4 (solid); (b) log-linear plot and (c) log-log plot of ^max vs. cq 
compared with the curves e^max/kT = VlOmM/co (dashed) and = 5 lnio(10mM/co) (soHd); and (d) log-log plot of ^max 
vs. Vmax compared to ^max = Vmax (solid). Points from the same experiment (varying concentration or voltage) are 
connected by line segments. 
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Figure 3: Typical experimental data (included in the estimates of Table 1) for two different types of nonlinear, induced- 
charge electrokinetic phenomena showing qualitative features not captured by the Standard Model, or the underlying 
electrokinetic equations of dilute solution theory, (a) Velocity of ACEO pumping of dilute aqueous solutions of KCl 
around a microfluidic loop by an asymmetric planar Au electrode array with the geometry of Refs 1 34 38 1 versus AC 
frequency at constant voltage, 3 Volts peak to peak (Vmax = 1.5 V), reproduced from Ref. |41 1. The data exhibits 
the unexplained flow reversal at high frequency (10 - 100 kHz) and strong concentration dependence first reported in 
Ref. |38|. (b) Velocity of ICEP motion of 5.7 //m metallo-dielectric Janus particles versus field-squared at diff'erent 
concentrations of NaCl in water at constant IkHz AC frequency, reproduced from Ref. |25 1. The data shows a similar 
decay of the velocity with increasing bulk salt concentration, which becomes difficult to observe experimentally above 
lOmM, in both experiments. 



2.2.2. Electrolyte dependence 

In addition to overestimating experimental velocities, the standard model fails to predict some 
important phenomena, even qualitatively. For example, ICEO flows have a strong sensitivity 
to solution composition, which is under-reported and unexplained. Most experimental work has 
focused on dilute electrolytes ll34ll3nil56[[2Ql . (See Table[T]) Some recent experiments suggest a 
logarithmic decay of the induced electro-osmotic mobility, b oc ln(Cc/co), with bulk concentration 
Co seen in KCl for ACEO micropumps |[38||4B|, in KCl and CaCOs for ICEO flows around metal 
posts 1 157 1, and in NaCl for ICEP motion of metallo-dielectric Janus particles | 25 1. This trend is 
visible to some extent at moderate concentrations in Fig.[2jb) over a wide range of experimental 
conditions, although a power-law decay also gives a reasonable fit at high salt concentrations. 
Two examples of diff'erent nonlinear electrokinetic phenomena (ACEO fluid pumping and ICEP 
particle motion) showing this unexplained decay with concentration are shown in Fig. [3] 

In all experiments, such as those in Fig. [2j the flow practically vanishes for cq > 10 mM, 
which is two orders of magnitude below the salinity of most biological fluids and bufl'er solu- 
tions (Co > 0.1 M). Experiments with DC |T58l|T59l and AC |[l60l[I^ field-eff'ect flow control, 
where a gate voltage controls the zeta potential of a dielectric channel surface, have likewise 
been limited to low salt concentrations below 10 mM in a variety of aqueous solutions. Sodium 
carbonate/bicarbonate bufl'er solutions have also been used in ICEO mixers with platinum struc- 
tures 1 141 1, but experimental data was only reported at a low ionic strength of 2.5 mM after 
dilution by water, and not for the original 1 M solutions. Remarkably, out of all the experimental 
work reviewed in this section, we could not find any observations of induced-charge electroki- 
netic phenomena at salt concentrations above 30 mM in water. 

12 



The Standard Model seems unable to explain the decay of flow with increasing salt concen- 
tration quantitatively, although it does aid in qualitative understanding. Substituting the Debye- 
Hiickel screening length for a binary z : z electrolyte in ^ we obtain 



1 



1 + V^oTq 

where 



for Co » Cc (8) 



kT I ss \ SbkT 



2sb \hszej 2(ze)^Aj 

is a crossover concentration, above which the flow decays like the inverse square-root of concen- 
tration. As noted above, it is common to attribute the theoretical over-prediction of ICEO flows, 
even in very dilute solutions to a large voltage drop across the compact layer (S » 1), but this 
would imply a strong concentration dependence (cq » Cc) that is not observed. Alternatively, 
fitting the compact-layer capacitance to reproduce the transition from dilute to concentrated so- 
lution behavior (cq ~ c^, ^ ~ 1) would eliminate the correction factor in dilute solutions (S ^ 1), 
making the theory again over-predict the observed velocities. For example, such difficulties are 
apparent in Ref. 1551 where this argument applied to the data of Gangwal et al 1251 for ICEP 
motion of metallo-dielectric Janus particles (Fig.|3ja)). 

Beyond the dependence on salt concentration, another failing of dilute-solution theory is the 
inability to explain the experimentally observed ion- specificity of ICEO phenomena. At the 
same bulk concentration, it has been reported that ICEO flow around metal posts [|157J , ACEO 



pumping by electrode arrays |41 1 and AC-field induced interactions in colloids 11611 depend on 
the ions. Comparing experiments under similar conditions with diff'erent electrolytes or diff'erent 
metal surfaces further suggests a strong sensitivity to the chemical composition of the double 
layer, although more systematic study is needed. In any case, none of these eff'ects can be cap- 
tured by the Standard Model, which posits that the ions are simply mathematical points in a 
dielectric continuum and that the surface is a homogeneous conductor or dielectric; all specific 
physical or chemical properties of the ions, solvent molecules, and the surface are neglected. 

2.2.3. Flow reversal 

In many situations of large induced voltages, the Standard Model does not even correctly 
predict the direction of the flow, let alone its magnitude. Flow reversal was first reported around 
tin particles in water |30|, where the velocity agreed with the theory |[T8l|28l sketched in Fig. 
la only for micron-sized particles and reversed for larger ones (90 - 400/im). The transition 
occured when several volts was applied across the particle and reversal was conjectured to be 
due to Faradaic reactions |30|. In this regime, reverse flows have recently been observed around 
large (millimeter scale) copper washers and steel beads with flow patterns resembling second- 
kind electro-osmosis II121II (see below); although the field was kept below the level causing 
gas bubbles at the anodic side of the metal, copper dendrites (resulting from electrodeposition) 
were observed on the cathodic side in dilute CuCl2 solutions, implying normal currents and 
concentration gradients. 

In microfluidic systems, flow reversal has also been observed at high voltage (> 2 V) and high 
frequency (10-100 kHz) in ACEO pumping by lOyum-scale planar electrode arrays for dilute 
KCl (381 Sni El (SU, as shown in Fig. [3jb), although not for water in the same pump geome- 
try f34j |40l. Non-planar 3D stepped electrodes |44 | can be designed that do not exhibit flow 
reversal, as demonstrated for KCl 1.45 J and water I48 L but certain non-optimal 3D geometries 
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can still reverse, as shown in water f40l . In the latter case the frequency spectrum also develops 
a double peak with the onset of flow reversal around 3 Volts peak to peak. In travelling-wave 
electro-osmosis (TWEO) in aqueous electrolytes ||49l|50|, strong flow reversal at high voltage 
has also been observed, spanning all frequencies |[50ll5n[T62 L and not yet fully understood. 

Flow reversal of ACEO was first attributed to Faradaic reactions under diff'erent condi- 
tions of larger voltages (8-14 V) and frequencies (1-14 MHz) in concentrated NaCl solutions 
(0.001 - 0.1 S/m) with a lOOyum-scale T-shaped electode pair composed of diff'erent metals (Pt, 
Al, Chromel) [163 1. Indeed, clear signs of Faradaic reactions (e.g. gas bubbles from electrolysis 
of water) can always be observed at sufficiently large voltage, low frequency and high concen- 
tration (381 11631 . In recent TWEO experiments |162, 164], signatures of Faradaic reactions 
(including pH gradients from electrolysis) have been correlated with low-frequency flow rever- 
sal at high voltage and bulk electroconvection has been implicated 11641 (see below). Under 
similar conditions another possible source of flow reversal is AC electrothermal flow driven by 
bulk Joule heating 1 165 L which has been implicated in reverse pumping over planar electrode 
arrays at high salt concentrations 11661 . Closer to standard ACEO conditions, e.g. at 1-2 V and 
50-100 Hz in water with Au electrode arrays, flow reversal can also be induced by applying a 
DC bias voltage of the same magnitude as the AC voltage L 1671 1 1681 1 1691 . Reverse ACEO flow 
due to "Faradaic charging" (as opposed to the standard case of "capacitive charging") is hypoth- 
esized to grow exponentially with voltage above a threshold for a given the electrolyte/metal 
interface 1 163, 168], but no quantitative theory has been developed. 

Simulations of the standard low-voltage model with Butler- Volmer kinetics for Faradaic reac- 
tions have only managed to predict weak flow reversal at low frequency in ACEO (121 Il2l ESI 
and TWEO [52^, 151 1 . In the case of ACEO with a planar, asymmetric electrode array, this effect 
has recently been observed using sensitive (yum/s) velocity measurements in dilute KCl with Pt 
electrodes at low voltage (< 1.5 V) and low frequency (< 20 kHz) |170|. Faradaic reactions can 
also produce an oscillating quasi-neutral diffusion layer between the charged diffuse layer and 
the uniform bulk, due to the normal flux of ions involved in reactions, and it has been shown 
via a low- voltage, linearized analysis of TWEO that flow reversal can arise in the case of ions 
of unequal diffusivities due to enhanced diffusion-layer forces on the fluid 1 151 1. Recently, flow 
reversal has been successfully predicted in TWEO by such a model allowing for bulk electro- 
convection in regions of pH gradients from electrolysis reactions and compared to experimental 
data |164|. However, the theory is still incomplete, and it seems that current models cannot 
predict the strong (> lOOyum/s), high-frequency (> 10 kHz) flow reversal seen in many ACEO 
and TWEO experiments | 38] |40l El [411. Faradaic reactions generally reduce the flow at low 
frequency by acting as a resistive pathway to "short circuit" the capacitive charging of the double 
layer (42|[52l, and diffusion-layer phenomena are also mostly limited to low frequency. Resolv- 
ing the apparent paradox of high-frequency flow reversal is a major motivation for our study. 

2.3. Nonlinear dynamics in a dilute solution 

Dilute- solution theories generally predict that nonlinear efl'ects dominate at low frequency. 
One reason is that the differential capacitance Co of the diffuse layer, and thus the "RC" time for 
capacitive charging of a metal surface, grows exponentially with voltage in nonlinear Poisson- 
Boltzmann (PB) theory. The familiar PB formula for the diffuse-layer differential capacitance of 
a symmetric binary electrolyte li6i .l29.L 
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2kT 



(10) 



was first derived by Chapman [171], based on Gouy's solution of the PB model for a flat diff'use 
layer 1 172|. It has been shown that this nonhnearity shifts the dominant flow to lower frequencies 
at high voltage in ACEO |42 1 or TWEO 1 173 1 pumping. It also tends to suppress the flow with 
increasing voltage at fixed frequency, since there is not adequate time for complete capacitive 
charging in a single AC period. 

At the same voltage where nonlinear capacitance becomes important, dilute- solution theory 
also predicts that salt adsorption |[T29l fT74l IT3T1 [TTSl fTTSl and tangential conduction [1311 [ml 
by the diff'use layer also occur and are coupled to (much slower) bulk diff'usion of neutral salt, 
which would enter again at low frequency in cases of AC forcing. If concentration gradients have 
time to develop, then they generally alter the electric field ("concentration polarization") and can 
drive bulk electroconvection |164|, diff'usio-osmotic slip L177. ,.8. .123J and in some cases, non- 
equilibrium space charge and second-kind electro-osmotic flow 11161 11181 11171 11241 11761 (if the 
bulk concentration goes to zero, at a limiting current). 

Bulk gradients in salt concentration may play a crucial role in induced-charge electrokinet- 
ics, especially at voltages large enough to drive Faradaic reactions. Concentration polarization 
has been demonstrated around electrically floating and (presumably) blocking metal posts in DC 
fields and applied to microfluidic demixing of electrolytes 1 1781. In nonlinear electrokinetic 
theory, difl'usion-layer phenomena have begun to be considered in low-voltage, linearized anal- 
ysis of TWEO with Faradaic reactions 1 151 1, and flow reversal of TWEO has been attributed 
to bulk electro-convection in regions of concentration polarization, for asymmetric electrolytes 
with unequal ionic mobilities 11641 . Such efl'ects could be greatly enhanced in the strongly non- 
linear regime and as yet unexplored. Including all of these efl'ects in models of induced-charge 
electrokinetic phenomena presents a formidable mathematical challenge. 

To our knowledge, such complete nonlinear modeling within the framework of dilute solution 
theory has only begun to be accomplished in the case of ACEO pumping (albeit without Faradaic 
reactions) in the Ph.D. thesis of Olesen [43 1 and the papers of Suh and Kang [1751 1179|[ by 
applying asymptotic boundary-layer methods to the classical electrokinetic equations in the thin- 
double-layer limit. At least in this representative case, all of the nonlinear large- voltage efl'ects 
in dilute solution theory in the solution phase tend to make the agreement with experiment worse 
than in the Standard Model [43 |. The flow is greatly reduced and shifts to low frequency, while 
the efl'ects of salt concentration and ion-specificity are not captured. Similar conclusions have 
been reached by a recent numerical and experimental study of fixed-potential ICEO for DC bias 
of 9 Volts [11461 . where the correction factor is found to be A = 0.005 for the linear theory, but 
only A = 0.05 if nonlinear capacitance (fTOl) and surface conduction from PB theory are included 
in the model (albeit without accounting for bulk concentration gradients). 

On the other hand, Suh and Kang [11791 have recently shown that the experiments of Green et 
al. on ACEO over a symmetric electrode pair (33] can be better described by nonlinear dilute 
solution theory (without surface conduction, but with oscillating diff'usion layers), if the no-flux 
boundary condition for ideally polarizable electrodes is replaced by a model of surface adsorption 
of ions [|175i using a Langmuir isotherm applied at the Stern plane fl80l ll8li . This approach re- 
lies on the classical, but somewhat arbitrary, partitioning of the double layer into an outer difl'use 
part (where "free" ions can drive fluid flow, governed by continuum equations of dilute solution 
theory) and an inner compact part (where "adsorbed" ions cannot, described by boundary con- 
ditions). It makes sense to describe ions that make contact with the electrode surface (breaking 
free of their own solvation shells and penetrating that of the metal, i.e. jumping from the "outer 
Helmholtz plane" to the "inner Helmholtz plane" (155 \) via an adsorption boundary condition, 
but the transition to dilute solution behavior is surely not so abrubt. Indeed, we will see that 
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concentrated solution theories of the Hquid phase can predict an effectively immobile compact 
layer forming at high voltage and advancing into the solution. 

Certainly, more theoretical work is needed on nonlinear dynamics of electrolytes in response to 
large voltages, especially in the presence of Faradaic reactions, but we believe the time has come 
to also question the validity of the underlying electrokinetic equations themselves. Based on 
experimental and theoretical results for induced-charge electrokinetic phenomena, we conclude 
that dilute solution theories may not properly describe the dynamics of electrolytes at large volt- 
ages. In the following sections, we consider some fundamental changes to the Standard Model 
and the underlying electrokinetic equations, while preserving the mean-field and local-density 
approximations, which permit a simple mathematical description in terms of partial diff'erential 
equations. We review relevant aspects of concentrated-solution theories and develop some new 
ideas as well. Through a variety of model problems in nonlinear electrokinetics, we make the- 
oretical predictions using modified electrokinetic equations, which illustrate qualitatively new 
phenomena, not predicted by the Standard Model and begin to resolve some of the experimental 
puzzles highlighted above. 



3. Crowding effects in a concentrated solution 

3.1. Mean-field local-density approximations 

3.1.1. Modified Poisson-Boltzmann theories 

All dilute solution theories, which describe point-like ions in a mean-field approximation, 
break down when the crowding of ions becomes significant, and steric repulsion and correlations 
potentially become important. If this can be translated into a characteristic length scale a for 
the distance between ions, then the validity of Poisson-Boltzmann theory is limited by a cutoff' 
concentration c^ax = ^"^^ which is reached at a fairly small diff'use-layer voltage, 

v}/^ ^ _^ ln(^L ^ Inia'co). (11) 
ze \ Co I ze 

where z is the valence (including its sign) and cq the bulk concentration of the counterions. In 
a dilute solution of small ions, this leads to cutoff's well below typical voltages for ICEO ffows. 
For example, even if only steric eff'ects are taken into account, with e.g. a = 3 A(for solvated 
bulk - Cr interactions 1 182|), then ~ 0.33V for cq = lO'^ M and z = 1. 

To account for the obvious excess ions in PB theory. Stern 1 1831 long ago postulated a static 
compact monolayer of solvated ions 1 155 |. A similar cutoff' is also invoked in models of ICEO 
ffows, where a constant capacitance is added to model the Stern layer and/or a dielectric coating. 



which carries most of the voltage when the diff'use-layer capacitance ([TO]) diverges. However, it 
seems unrealistic that a monolayer could withstand most of the voltage drop in induced-charge 
electrokinetic phenomena at a blocking surface (e.g. without dielectric breakdown [1841). In 
any case, a dynamical model is required for a "condensed layer" that is built and destroyed as the 
applied field alternates. As sketched in Fig.|4j the condensed layer forms in the diff'use part of the 
double layer and thus should be described by the same ion transport equations. For a non-ideally 
blocking surface, it may also be necessary to account for surface adsorption of ions (breaking 
free of solvation shells) and Faradaic reactions (electron transfer) via compact-layer boundary 
conditions (see below), but these should be applied to a model of the diff'use layer that allows for 
the crowding of ions near the surface, which must occur to some degree with increasing voltage. 
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Figure 4: Sketch of solvated counterions (larger green spheres) and co-ions (smaller orange spheres) near a polarizable 
surface, (a) At small induced voltages, ^ the neutral bulk is only slightly perturbed with a diffuse-charge layer 
of excess counterions at the scale of Ad- (b) At moderate voltages, ~ the diffuse layer contracts, as described by 
Poisson-Boltzmann (PB) theory, (c) At large voltages, ^ the counterions inevitably become crowded, causing 
expansion of the diffuse layer compared to the predictions of the classical Gouy-Chapman-Stern model, sketched in (d), 
which is based PB theory for point-like ions with a minimum distance of approach, the "outer Helmholtz plane" (OHP), 
to model solvation of the surface. 
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A plethora of "modified Poisson-Boltzmann" (MPB) theories have been proposed to describe 
equifibrium ion profiles near a charged wall. As described in recent reviews |2Q1 |62j 11851 11861 
11871 11881 11891 ), there are many possible modifications to describe diff'erent physical eff'ects, 
such as dielectric relaxation, electrostatic correlations and volume constraints. In this paper, 
we focus on the simplest continuum models, which are based on the local-density and mean- 
field approximations. In spite of various limitations discussed below, this class of models is 
very convenient (if not required) for mathematical analysis and numerical simulation of time- 
dependent nonlinear problems (based on modified Poisson-Nernst-Planck equations 111 301 [75l 

ESI). 

The starting point for continuum modeling is a theory for the excess electrochemical potential 
of an ion 

^4'=^ii-^if'^' = kT\nfi, (12) 

relative to its ideal value in a dilute solution, 

iJf'^^ = kT\nci+Zie(p, (13) 

where Ci is the mean concentration and / is the chemical activity coefficient. (Equivalently, one 
can write /// = kTln^Xi) + ztecp, where Ai = fiCi is the absolute chemical activity |190|.) In the 



mean-field approximation (MF), the electrostatic potential in ([13]) self-consistently solves the 
MPB equation, 

-V-(£V0)=p = ^Z/^c„ (14) 

i 

where the source of the electric field acting on an individual ion is the mean charge density p, 
rather than the sum of fluctuating discrete charges. Time-dependent modified PNP equations 



then express mass conservation with gradient-driven fluxes I13QL as described below. 

In the asymptotic limit of thin double layers, it is often justified to assume that the ions are in 
thermal equilibrium, if the normal current is not too large and the nearby bulk salt concentration 
is not too low I154[|19lll , even in the presence of electro-osmotic flow I123II124L In terms of 
electrochemical potentials, the algebraic system {yU/ = constant} then determines the ion profiles 
Ci in the diff'use layer, which lead to eff'ective surface conservation laws 1 114J . In dilute-solution 
theory {j/^ - 0), this procedure yields the Boltzmann distribution, 

c,#) = cOexp(^), (15) 

where = - 0/, is the potential relative to its bulk value 0/^ just outside the double layer. For a 



symmetric binary electrolyte (z± = ±z), substituting into ([14]) yields the standard form of the PB 
equation 

(^00' = -P(0) = 2coZ^ sinh (^) , (16) 

from the Gouy-Chapman model of the double layer 111721 11711 . where s is typically set to a 
constant bulk value of E}y (which we relax below). In that case, by scaling length to and the 
potential to kT jze, the PB equation takes the simple dimensionless form, xf/" - sinh0. 

In concentrated- solution theories, the simplest and most common approach is based on the 
local density approximation (LDA), where depends only on the local ion densities, in the 
same way as in a homogeneous bulk system l6^ l62l[T85lll92|[T93l . The choice of a model for 
yu^^ then yields a modified charge density profile p, diflfering from ( [T6] ) with increasing voltage. 
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In this section, we focus on entropic effects in bulk lattice-gas and hard- sphere models, where 
yu^^ depends only on the local ion concentrations, but not, e.g., explicitly on the distance to a 
wall I1Q61 11941 or non-local integrals of the ion concentrations 11951 11961 1197L as discussed 
below. 

For any MF-LDA model, the equilibrium charge density can be expressed as a function of 
the potential, p(</^), but the Boltzmann distribution ([16]) is modified for non-ideal behavior, as we 
now demonstrate. Here, we set the permittivity to its bulk value, s = =constant, but below 



we will extend the derivation for field-dependent dielectric response £(E) in section 3.1.6 Given 
p(</^), we integrate the MPB equation (multiplied by i//') to obtain the electrostatic pressure, Pe 
and normal electric field at the inner edge of the diff'use layer. 



Pe('i'D)=\sEl= f pW#. (17) 

From the total diff'use charge per unit area, 

q = -signCi^D)^2speCyDX (18) 
we then arrive at a general formula for the diff'erential capacitance. 



which reduces to ([TO]) in a dilute solution. 

3.1.2. The Bikerman-Freise formula 

Since most MPB models are not analytically tractable, we first illustrate the generic conse- 
quences of steric eff'ects using the oldest and simplest mean-field theory 11831 l59l 11981 [76l [6Ql 
[8T1l . This model has a long and colorful history of rediscovery in diff'erent communities and 
countries (pieced together here with the help of P. M. Biesheuvel, Wageningen). It is widely 
recognized that O. Stern in 1924 11831 was the first to cutoff' the unphysical divergences of the 
Gouy-Chapman model of the double layer [1721 11711 by introducing the concept of a "com- 
pact layer" or "inner layer" of solvent molecules (and possibly adsorbed ions) forming a thin 
mono-molecular coating separating an electrode from the "diff'use layer" in the electrolyte phase. 
The resulting two-part model of the double layer has since become ingrained in electrochem- 
istry 1 155 1 . Over the years, however, it has somehow been overlooked that in the same ground- 
breaking paper I183L Stern also considered volume constraints on ions in the electrolyte phase 
and in his last paragraph, remarkably, wrote down a modified charge- voltage relation [his Eq. 



(2')] very similar to Eq. ([24]) below, decades ahead of its time. We have managed to find only 
one reference to Stern's formula, in a footnote by Freise 11981 . 

Although Stern had clearly introduced the key concepts, it seems the first complete MPB 
model with steric eff'ects in the electrolyte phase was proposed by J. J. Bikerman in 1942, 
in a brilliant, but poorly known paper |[59]| . (Bikerman also considered polarization forces 
on hydrated-ion dipoles in non-uniform fields, discussed below.) Over the past sixty years, 
Bikerman's MPB equation has been independently reformulated by many authors around the 
world, including Grimley and Mott (1947) in England EH |86l, Dutta and Bagchi (1950) in 
India |[l99l [200l EqU [SI, Wicke and Eigen (1951) in Germany ||203l [2041 |20l, Strat- 
ing and Wiegel (1993) in the Netherlands EoSl |207l [H i2Q8J, IgUc and Kralj-Ighc (1994) 
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in Slovenia 12091 f2T0', "ITT, "2T2l, and Bomkhov, Andelman and Orland (1997) in Israel and 
France (761 12131 [2141 . For an early review of electrolyte theory, which cites papers of Dutta, 
Bagchi, Wicke and Eigen up to 1954 (but not Bikerman or Freise, discussed below), see Redlich 
and Jones 12151 . 

Unlike Bikerman, who applied continuum volume constraints to PB theory, most of these 
subsequent authors derived the same MPB model starting from the bulk statistical mechanics 
of ions and solvent molecules on a cubic lattice of spacing a in the continuum limit (MF- 
LDA) where the concentration profiles vary slowly over the lattice. (More recently, "semi- 
discrete" lattice-gas models have also been developed, which consist of discrete layers de- 
scribed via mean-field approximations, without taking the continuum limit in the normal direc- 
tion [11881 [87112161 1217II .) While early authors were concerned with departures from PB theory 
in concentrated electrolytes 1591 1 199(1205 1 or ionic crystals |[85][86l, recent interest in the very 
same mean-field model proposed by Bikerman has been motivated by a wide range of modem 
applications involving electrolytes with large ions or biological molecules | 74l [76l 12131 1214L 
polyelectolytes 12181 12191 1220L polymeric electrolytes I1221L electromagnetic waves in elec- 
trolytes |75 |, electro-osmosis in nanopores |68 , 69], electrophoresis of colloids |[7Ql[7n[72ll73l , 
and solvent- free ionic liquids (81] [82] [8^, in addition to our own work on dilute electrolytes in 
large applied voltages (50| [T30| [58l . 

In the present terminology, Bikerman's model corresponds to an excess chemical potential 



IJi^^ = -kT ln(l - O) (Bikerman) , (20) 

associated with the entropy of the solvent, where O = YjI Ci is the local volume fraction of 
solvated ions on the lattice (62l . For now, we also assume a symmetric binary electrolyte, c+ = 
= Co, z± = ±z, to obtain an analytically tractable model. As shown in Fig.|5[a), when a large 
voltage is applied, the counterion concentration exhibits a smooth transition from an outer PB 
profile to a condensed layer at c = Cmax = cT^ near the surface. Due to the underlying lattice-gas 
model for excluded volume, the ion profiles eff'ectively obey Fermi-Dirac statistics, 

c± = 5 , (21) 

1 +2v^m\\^{zeil/l2kT) 

where v = 2a^co = ^buik is the bulk volume fraction of solvated ions. Classical Boltzmann 
statistics and the Gouy-Chapman PB model are recovered in the limit of point-like ions, v = 0. 

For a flat double layer, similar results can be obtained with the even simpler Composite Difl'use 
Layer model of Kilic et al. 16^ (also termed the "cutofl' model" in Ref. 11851 ), where an outer 
PB difl'use layer is abrubtly patched with an inner condensed layer of only counterions at the 
uniform, maximal charge density. This appealingly simple construction requires assumptions 
about the shape of the condensed layer (e.g. a plane), so its position can be determined only 
from its thickness or total charge. Even if it can be uniquely defined, the cutoff' model introduces 
discontinuities in the co-ion concentration (which drops to zero in the condensed layer) and in 
the gradient of the counter-ion concentration, although the same is also true of Stern's original 
model of the compact layer. In this work we focus on Bikerman's model since it is the simplest 
general model of steric eff'ects that remains analytically tractable; unlike the cutoff' model, it 
predicts smooth ionic concentration profiles in any geometry and can be naturally extended to 
time-dependent problems [130! 751 . 

Substituting the equilibrium ion distributions ( [2T] ) into Poisson's equation ([14]), we obtain 
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Figure 5: (a) The equilibrium distribution of counterions in a flat difl'use layer for large applied voltages zeWo/kT = 
10, 20, ... , 100 predicted by Poisson-Boltzmann theory (PB) and Bikerman's modified theory (MPB) taking into account 
an eff'ective (solvated) ion size a, where v = la^co = 0.001 is the bulk volume fraction of solvated ions, (b) The 
diff'use-layer diff'erential capacitance C/) vs. voltage predicted by PB fTo) (y = 0) and MPB {25) (y > 0), scaled to the 
low- voltage Debye-Hiickel limit s/Ad(co). 



Bikerman's MPB equation, 



2coze smh(zei// / kT) 



(22) 



1 +2ysinh^(z^j/^/2^r)' 
which has the simple dimensionless form (scahng length to Ad, potential to kT/ze), 

~„ ^ sinh<A 

1 +2ysinh2(jA/2)' 

extending the PB equation for a nonzero volume fraction y of ions in the bulk solution. Integrat- 



(23) 



ing across the diffuse layer, the charge- voltage relation (18) takes the form. 



qy = sgn(^z))2z^<^o^D J - In 



1 + 2y sinh^ 



2kT 



(24) 



which was probably first derived by Grimley L86J in a lattice-gas theory of diff'use charge in 
ionic crystals 1851 , independent of Bikerman. Grimley 's formula has the same form as Stern's 
surprising Eq. (2') noted above 1 183|, but it has all the constants correct and clearly derived. 
Recently, Soestbergen et al.| 221 1 have given a slightly diff'erent formula approximating ([24]) that 
is easier to evaluate numerically for large voltages, and applied it to ion transport in epoxy resins 
encapsulating integrated circuits. 

Although Stem and Grimley derived the modified form of the charge-voltage relation with vol- 
ume constraints, they did not point out its striking qualitative diff'erences with Gouy-Chapman 
dilute- solution theory, noticed by several recent authors |[60||8T1. This important aspect was ap- 



parently first clarified by Freise 1 198 L who took the derivative of (24 ) and derived the diff'erential 



capacitance ( 19 ) in the form 



^ r 



sinh( 



kT ) 



[1 + 2vsinh2 (f^)] ^^ln[l+2vsinh2(f^)l 
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(25) 



and pointed out that decays at large voltages. This elegant formula, also derived by Kilic 
et al. 1 60 1 and Kornyshev | 81 1 (and Oldham | 84| for the ionic-liquid limit v = 1), predicts the 



opposite dependence of Chapman's formula (10) from dilute- solution theory, which diverges ex- 



ponentially with \^d\- Since Chapman 1 171 1 is given credit in "Gouy-Chapman theory" for first 



deriving the capacitance formula ([TO]) for Gouy's original PB model 1 172], we suggest calling 
Eq. (25) the "Bikerman-Freise formula" (BF), in honor of Bikerman, who first postulated the 
underlying MPB theory, and Freise, who first derived and interpreted the modified diff'erential 
capacitance. By this argument, it would be reasonable to also refer to the general MPB model 
as "Bikerman-Freise theory", but we will simply call it "Bikerman's model" below, following 
Refs. |[68l[69l|6l. 



According to the BF formula (25 ), as shown in Fig.|5jb), the differential capacitance increases 



following PB theory up to a maximum near the critical voltage ^ (^^d then decreases as 
the square-root of the voltage, 

because the eff'ective capacitor width grows due to steric eff'ects, as seen in Fig. [JJa). In stark 
contrast, the PB diff'use-layer capacitance diverges exponentially according to Eq. ([TO]), since 
point-like ions pile up at the surface. Although other eff'ects, such as compact layer compression 
and specific adsorption of ions 1 1881 1155 1 (discussed below) can cause the diff'erential capac- 
itance to increase at intermediate voltages, this eff'ect is quite general. As long as the surface 
continues to block Faradaic current, then the existence of steric volume constraints for ions im- 
plies the growth of an extended condensed layer at sufl&ciently large voltages, and a concomitant, 
universal decay of the diff'erential capacitance. Indeed, this eff'ect can be observed for interfaces 
with little specific adsorption, such as NaF and KPFe on Ag 1222112231 or Au 12241 [2251 [2261 . 
(See Fig. [8] below for fits to simple mean-field models.) The same square-root dependence at 
large voltage can also be observed in experiments | 81 1 and simulations |[82j[83l of ionic liquids 
at blocking electrodes, with remarkable accuracy. We conclude that the decay of the double-layer 
diff'erential capacitance at large voltage is a universal consequence of the crowding of finite-sized, 
mobile charge carriers near a highly charged, blocking surface. 

The BF formula ([26]) also illustrates another general feature of double-layer models with steric 
constraints, shown in Fig. [6] The differential capacitance at large voltages is independent of bulk 
salt concentration, hut ion specific through z and a. This prediction is reminiscent of Stern's 
picture picture 1 183 J of an inner, compact layer of solvent molecules which carries the majority of 
a large double layer voltage, compared the outer, diff'use layer described by dilute PB theory. This 
picture is ubiquitous in electrochemistry 1 155 1, and first gained acceptance based on Grahame's 
famous experiments on mercury drop electrodes 12271 12281 [T88J . The significant diff'erence, 
however, is that the condensed layer forms continuously in the solution near the inner edge of 
the diff'use layer due to ion crowding eff'ects in a general model of the electrolyte phase, which is 
not restricted to flat quasi-equilibrium double layers. 



3.1.3. Hard- sphere liquid models 

Although Bikerman's model describes steric efl'ects in a convenient and robust analytical form, 
the bulk ionic volume fraction y is best viewed as an empirical fitting parameter. For crystalline 
solid electrolytes, its microscopic basis in a lattice model is realistic, but even then, the thinness 
of the condensed layer, comparable to the lattice spacing at normal voltages, calls the continuum 
limit into question. For liquid electrolytes involved in electrokinetic phenomena, it would seem 
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Figure 6: Differential capacitance C/) vs. voltage in Bikerman's MPB model {25) with a = 4 Afor y values corresponding 
to Co = 1, 10, 100 mM. In contrast to Fig.[5jb), here Co is scaled to a single constant, ImM), for all concentrations. 



more realistic to start with the "restricted primitive model" of charged hard spheres in a uniform 
dielectric continuum 12291 in developing better MPB models Il6nil92[|62ll . From this theoreti- 
cal perspective, Bikerman's lattice-based model has the problem that it grossly under-estimates 
steric effects in hard-sphere liquids; for example, in the case of a monodisperse hard- sphere 
liquid, the volume excluded by a particle is eight times its own volume I23QI 12311 [62l . Al- 
though we focus on electrolytes at large voltages, it is also interesting to consider the MF-LDA 
dynamics of charged hard spheres to model other systems, such as dense colloids 12321 1233L 
poly electrolytes 12341 [235ll . and ionic liquids l8Tll82ll. 

Various MF-LDA approximations of ji^.^ for bulk hard- sphere liquids can be used to develop 
more sophisticated MPB models, which yield similar qualitative behavior of the diffuse-layer 
differential capacitance 11921 [621 , due to the generic arguments given above. For example, the 
Carnahan- Starling (CS) equation of state for a bulk monodisperse hard- sphere liquid corresponds 
to the following excess chemical potential 12361 [Til, 

T¥ = ^^^^3 (Carnahan-Starling) (27) 

kT (1 - (^y 

Although this algebraic form precludes analytical results, it is much simpler to evaluate nu- 
merically and incorporate into continuum models of electrokinetic phenomena than are more 
sophisticated MPB approximations which go beyond the LDA, e.g. based on self-consistent cor- 
relation functions |78, 79, 80, 237 1 or density functional theory |238, 239, 240 1, which require 
solving nonlinear integro-differential equations, even for a flat double layer in equilibrium. As 
shown in Fig.|7jb), the simple CS MPB model predicts capacitance curves similar to Fig. [6] with 
Bikerman's model, respectively, only with more realistic salt concentrations |62|. In particular, 
the differential capacitance in Bikerman's model ressembles that of CS MPB if an unrealistically 
large hydrated ion size a (or large bulk volume fraction y) is used, due to the under-estimation of 
liquid steric effects noted above. 

23 



(c) 




Figure 7: Modified Poisson-Boltzmann theory for a binary solution of charged hard spheres of diameter a = 4 Ausing 
the Carnahan-StarUng (CS) equation of state {27) . (a,b) The diffuse-layer differential capacitance vs. voltage, analagous 
to Fig. [5]b and Fig. [6] respectively, (c) The counterion density profile in the diffuse layer at voltages ze^o/kT = 
5, 10, 20, 40, 60, 80, 100 and concentration of Co = 10 mM, analogous to Fig.|5ja). (d) The surface counterion density vs. 
voltage at Co = 10 mM in the CS and Bikerman MPB models. 
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In spite of similar-looking capacitance curves, however, there are important differences in the 
ionic profiles predicted by the two models. As shown in Fig [Sja), in Bikerman's model steric 
effects are very weak until the voltage becomes large enough to form a thin condensed layer at 
maximum packing. As such, the width of the diffuse layer at typical large voltages is still an order 
of magnitude smaller than the Debye length Ad relevant for small voltages. In contrast, steric 
effects in a hard-sphere liquid are stronger and cause the diffuse layer to expand with voltage 
as shown in Fig. [TJc). The widening of the diffuse layer reduces its differential capacitance, 
but without forming the clearly separated condensed layer predicted by Bikerman's model. As 
shown in Fig. [TJd), the counterion density at the surface in the CS MPB model increases more 
slowly with voltage as compared to the Bikerman model. These differences will be important 
when we discuss the viscosity effects in Section 4. 

An advantage of the hard-sphere approach to volume constraints is that it has a simple exten- 
sion to mixtures of unequal particle sizes 112411 which can be applied to general multicomponent 
electrolytes |61, 62, 1921 12351 12341 . According to the Boublik-Mansoori-Camahan-Starling- 
Leland (BMCSL) equation of state 112421 12431 , the excess chemical potential of species / in a 
mixture of N species of hard spheres with different diameters {a/} is given by 

- ' 1 -h ^ ln(l-0)-h ' 



kT \ 03 (D2 j ' l-O 



(l-0)4§"H-^^^ U^(l-a))3 ^^^C^L) (28) 



where ^„ = ^JLi ^flY^ ^ is the volume fraction of species 7, and O = is the total 

volume fraction of ions. Although this formula may seem complicated, it is an algebraic ex- 
pression that can be easily expanded or evaluated numerically and thus is much simpler than 
statistical theories based on integral equations 112411 11861 1187II . The first BMCSL correction to 
dilute solution theory is simply, 

,,ex N I v3 

The BMCSL MPB model for asymmetric electrolytes predicts the segregation of ions of different 
size and/or charge in the diffuse layer 1621 and has been applied to adsorption phenomena in 
poly electrolyte layers Ii2351 12341 . The broken symmetry between ions of different sizes is an 
important qualitative effect, which we will show implies new electrokinetic phenomena at large 
voltages, regardless of the model. 

A word of caution: In spite of its mathematical convenience, the local-density approximation 
is known to provide a poor description of confined hard-sphere liquids, even in equilibrium. For 
example, it cannot capture density oscillations near a wall or two-point correlation functions of 
hard spheres. These features can be approximated by various weighted density approximations 
(WDA) in statistical mechanics, which redefine the local reference densities Ci as averages over 
the inhomogeneous densities with a suitable weight function 12441 12451 124^1 [2771 [2481 fT89^ 
The non-local weighted densities are used in place of the local densities in chemical potential 



expressions for homogeneous systems, such as (27) and (28) for hard spheres, and several pre- 



scriptions for the weight function are available. (See section 5.2 ) Non-local MF-WDA models 
have been applied to electrolytes to quantify excluded volume effects involving ions and sol- 
vent molecules 12491 [T95.. 196.. 197.1 . Recent Monte-Carlo simulations of hard-sphere counterion 
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Li+ 


1.20 


7.64 


4.2 


Na+ 


1.90 


7.16 


4.0 




2.66 


6.62 


3.8 


cr 


3.62 


6.64 


3.6 




Table 2: Comparison of the bare ion diameter in a crystalline solid, dx, with the effective solvated diameter ds in water 
from bulk transport measurements (akin to a Stokes diameter) 1 252 1 and the "hard-sphere" diameter (akin to a collision 
cross section) inferred from viscosity data d^ in dilute aqueous solution 1 253 1 for some common ions used in nonlinear 
electrokinetic experiments. The figure depicts an ion with its eff'ective hard-sphere and solvation shells, in red and blue 
respectively. Note that the eff'ective sizes ds and d^ in solution are much larger than the bare ion size dx and exhibit 
diff'erent trends. In the text, we argue that the appropriate eff'ective ion size a in our models of highly charged double 
layers may be approximated by ds, and possibly larger. 



profiles around a hard-sphere macro-ion have shown that LDA can perform even worse than 
dilute- solution PB theory, while WDA theories are able to fit the simulations well |250|. Even 
the simplest WDA, however, is a non-local continuum theory and thus requires solving non- 
linear integral equations for the equilibrium densities, and integro-partial differential equations 
for time evolution. Clearly, LDA-based partial diff'erential equations are much better suited for 
mathematical modeling, if they can capture enough of the essential physics in a given situation. 

3.1.4. Interpretation of the effective ion size 

In order to apply our modified electrokinetic equations, we stress that the effective diameter 
of a solvated ion in various MF-LDA theories is not simply related to its bare atomic size and 
can exhibit very different trends. Smaller bare ions tend to be more heavily solvated and there- 
fore have larger effective diameters |251|. Effective solvated ion sizes depend on the size and 
charge of the ions, the nature of the solvent, the ion concentration, and temperature - as well 
as the mathematical models used in their definitions. Table [2] compares bare ionic diameters 
in crystalline solids to effective solvated-ion diameters inferred from bulk properties [252] and 
"hard- sphere" diameters inferred from viscosity measurements 12531 . both in aqueous solutions. 
In these models used to interpret experimental data, the hard sphere radius is essentially a col- 
lision size, whereas the the effective solvated radius is an effective size for transport properties, 
similar to a Stokes radius. The effective solvated radius is generally larger than the hard sphere 
value. Both are greatly exceed the bare diameter and exhibit roughly opposite trends with the 
chemical identity of the ion. 

What is the appropriate effective ion size al Unlike the models used to infer the various ion 
diameters in Table |2j our models seek to capture crowding effects in a highly charged double 
layer, rather than in a neutral bulk solution (albeit still using MF-LDA theories). As such, it is im- 
portant to think of crowded counterions of the same sign and not a neutral mixture of oppositely 
charged ions (where our models reduce to the Standard Model in typical situations with dilute 
electrolytes). Below, we will argue that the crowding of like-charged counterions in large elec- 
tric fields leads to some different physical effects, not evident in the neutral bulk liquid. Among 
them, we can already begin to discuss solvation effects. In the bulk, ions cannot reach very high 
concentrations due to solubility limits, but a condensed layer of counterions cannot recombine 
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Figure 8: Fitting of simple MF-LDA models of finite-sized ions to the experimental data of Valette (Fig. 3 of Ref. 12221 ) 
for the diff'erential capacitance of the aqueous KPFg | Ag double layer. For all curves, the electrode area is scaled by 
a surface roughness factor R = 1.1, and the theoretical diff'use-layer capacitance Co from fT9) is taken in series with a 
constant Stern-layer capacitance Cs = 125//F/cm^. The only other adjustable parameter in each pane l is the eff'ective 
solvated ion size a. In (a), we use a hard-sphere radius of a = 4 Ain the Carnahan-Starling model \2l\ , while in (b) we 
use a lattice spacing ofa= 11 Ain Bikerman's MPB model |2o)-|25). 



and is unaffected by solubility (except for the possibility of electron transfer reactions near the 
surface). Moreover, like-charged ions cannot easily "share" a solvation shell and become com- 
pressed to the hard-sphere limit, since the outer surfaces of the polarized solvation shells have 
the same sign and yield electrostatic repulsion. Therefore, we propose that the "ion size" in our 
models is an effective solvated ion size at high charge density, which is much larger than the 
bare crystalline and hard-sphere ion sizes in Table [2] and may also exceed the solvated ion size 
inferred from bulk transport models. This physical intuition is borne out by the comparisons 
between theory and experiment below for nonlinear electrokinetics, although we will not claim 
to reach any quantitive molecular-level conclusions. 

3.1.5. Comparison with experiments on blocking surfaces without ion adsorption 

To illustrate these points, in Fig. [8] we fit our simple mean-field models to experimental data of 
Valette 1 222 1 for the differential capacitance of a standard electrochemical interface with negligi- 
ble surface adsorption of ions: aqueous KPF6 solution with single-crystal Ag electrodes. These 
experiments and others were cited by Di Caprio et al. |192[|193| as motivation to develop mean- 
field hard-sphere models for double-layer capacitance, but they concluded that "a quantitative 
comparison with the experimental results is not possible at this stage". Here, we proceed naively 
to directly test the predictions of two simple theories. To model the liquid electrolyte, we use 
the simplest CS and Bikerman MPB models, each of which has only one adjustable parameter, 
the effective ion size a. To model the surface (and admittedly, to aid in the fitting), we introduce 
two standard fitting parameters: (i) a "surface roughness factor" Rs multiplying the nominal 
electrode area, and (ii) a Stern-layer capacitance Cs in series with the diffuse-layer capacitance, 
ostensibly to model the dipolar electrode solvation layer, since we already take into account the 
crowding of finite-sized solvated ions. Unlike Valette II222I . who adjusts Rs for each measure- 
ment, we choose a fixed value of Rs = 1.1 for all the data, since this parameter should reflect 
fixed surface profile variations at scales larger than the inner part of the double layer (where 
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ion crowding occurs) and thus should not depend on concentration or voltage. Perhaps more 
seriously, we assume a constant Stern capacitance Cs for simplicity, in order to highlight what 
nonlinear trends can be captured by the MPB liquid model alone. This is rather different from the 
standard fitting procedure in electrochemistry dating back to Grahame 12271 II 55L which essen- 
tially defines the compact-layer capacitance empirically to account for all eff'ects not predicted 
by the dilute- solution PB model of the diff'use layer. 

First, we consider the CS excess chemical potential ( [27) ) for a hard-sphere liquid in the general 
MPB formula for the diff'erential capacitance ([19]). As shown in Fig. [8] a reasonably good fit is 
obtained using Cs = 125yuF/cm^ and an eff'ective ion diameter of (2 = 4 A, comparable to the 
various ion sizes (in particular for K^) inferred from completely diff'erent measurements. Near 
the point of zero charge and in dilute solutions (< 100 mM), the modified theory reduces to the 
Gouy-Chapman formula for a dilute solution ( [TQ| , and we obtain an excellent fit of the data, 
consistent with volumes of prior work in electrochemistry II 2271 12281 11881 11551 . At the highest 
salt concentration (cq =100 mM), the theory underestimates the capacitance, although this is 
could be improved by adjusting Cs for each salt concentration, as is commonly done. Of course, 
it is more interesting for us to focus on the regime of large voltage, where ion crowding occurs 
in the model. From Fig. [8ja) it is clear that the CS model cannot fit the fairly sudden decay 
of Cd beyond the maximum, although it predicts the qualitative loss of bulk salt concentration 
dependence (due to the dominance of near-surface crowding eff'ects, as described above). This is 
due to the gradual onset of crowding eff'ects for hard spheres in large fields, illustrated earlier in 
Fig. 13 

Next, we fit Bikerman's MPB model to the same experimental data by fixing the parameters 
of the surface, Rs = 1.1 and Cs = 125yuF/cm^, and varying only the ion-size parameter a. We 
might expect a worse fit because the model (based on a lattice gas) has less statistical justification 
in the liquid state, but the results in Fig[8jb) show a significant improvement in fitting the voltage 
dependence, due to the more sudden onset of steric eff'ects with increasing voltage versus the 
CS model (Fig. [7]). Overall, Bikerman's model provides an impressive fit to the concentration 
and voltage dependence of Co across the entire experimental range, considering that it has only 
one adjustable parameter a to describe the liquid (and we added two more, Rs and Cs , for the 
surface). 

On the other hand, the excellent fit by Bikerman's model in Fig. [8jb) is achieved with a sur- 
prisingly large value for the eff'ective ion size, a = lA nm, at least twice the typical solvated 
ion sizes in the bulk liquid from Table [2j Interestingly, very similar large values of a result from 
fitting Bikernan's model to completely diff'erent experiments on ACEO pumping by electrode 
arrays [61], as discussed below. These results may point to correlation eff'ects at high charge 
density or other neglected eff'ects (see below), which are fortuitously fitted well by Bikerman's 
model with an enlarged ion size. 

It would be interesting to extend our capacitance analysis to electrolytes with little ion adsorp- 
tion on Au surfaces | 224, 225, 226 1, since gold surfaces are widely used in experimental studies 
of induced-charge electrokinetic phenomena (including most of the entries Table 1). Given the 
successful fit in Fig. [8jb), it would also be interesting to perform ICEO experiments using sil- 
ver surfaces and KPFg solutions or other electrolytes with similar capacitance behavior, showing 
little surface adsorption of ions. By taking the diff'erential capacitance from fitting independent 
measurements, one can directly test other assumptions in the model, such as the HS slip formula 
(or its modifications, derived below). 
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3.1.6. Dielectric response in a concentrated solution 

We close our discussion of mean-field models of ion crowding by briefly considering nonlinear 
dielectric response. We have already mentioned that Bikerman, in his pioneering paper on efl'ects 
of finite ion size 1591 , also modeled the polarization of solvated ions by the local electric field, 
in terms of induced dipoles. A related eff'ect is the polarization of a Stern monolayer of solvent 
dipoles on a metal electrode, which plays a major role in classical models of the compact inner 
layer 12281 [T88 1. This approach was extended by Macdonald and Kenkel 1871 12161 [2171 to also 
describe the diff'use part of the double layer by postulating discrete layers of finite-sized dipoles 
separating layers of ions, where each layer is treated by a mean-field approximation. 

Closer in spirit to the continuous MF-LDA theories of this section, Abrashkin, Andelman and 
Orland |254| recently derived a "Modified Dipolar PB" equation for a symmetric binary elec- 
trolyte by considering an equilibrium lattice gas of both ions and dipoles of (the same) finite size, 
thus generalizing Bikerman's model (or its many reincarnations cited above, including Borukhov, 
Andelman, Orland |76 1) in a natural way to allow for a variable dipole density. In our notation, 
their dipolar MPB equation takes the form 

,„ . CdPo (0(poi/^'/kT)\ 2cozesmh(zei///kT) 
'"^ ^ ^ i D I = D 

where is the vacuum permittivity, Cd is the bulk density of dipoles of moment po on a lattice 
of spacing a given by Cd + 2co = a~^, where cq is the bulk salt concentration; the second term 
describes solvent polarization as a generalization of the classical expression for a dilute solution 
of dipoles (below), where the function 

_ ~, sinh£ ^ ~, coshE" sinh£ 

g(E) = ^j:(E) = ^ W- ^''^ 

relates the nonlinear polarization to the Langevin function £(x) = coth(x) - 1 / x, which arises 
in the limiting probability distribution of a sum of independent, randomly oriented unit vectors 
(also known as Rayleigh's random walk 12551 ): the third term is the free charge density, where 



the familiar PB expression (16) appears in the numerator; the eff'ect of volume constraints enters 



through a generalized "Bikerman factor", 

D = 2cocosh(p^]^cd'^^^^^ (32) 
\kT) Poi/r'/kT ^ ^ 

which rescales both the solvent polarization and the free charge density. In the limit of linear 



dielectric response, poi//' /kT 0, the dipolar MPB equation Eq. (30) reduces to Bikerman's 



MPB equation (22), where the permittivity has the eff'ective bulk value 



£b = so+irr7p- (33) 
3kT 



In the limit of point-like ions and dipoles a ^ 0, Eq. ( |30| ) reduces to the PB equation with the 
classical field-dependent permittivity, 

s(E) = so^—&( — ) (34) 

It would be interesting to develop this model further, test it against experiments, and eventually 
apply it to electrokinetics. In its present form, however, the dipolar MPB theory defies the 
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common wisdom in electrochemistry by predicting a large increase in permittivity in the inner 
part of the double layer, due to a rise in the dipole density near the surface [254] . Perhaps 
this discrepancy is due to the breakdown of the LDA I25QI already noted above, and different 
predictions might result from non-local WDA theories with dielectric relaxation 11951 12491 11961 
fT97l . 

Based on half a century of fitting experimental capacitance data (with admittedly simple mod- 
els), electrochemists have come to believe that that the local dielectric constant of an electrolyte 
is generally reduced in the inner part of the double layer, due to the alignment of solvent dipoles 
in the large local field ("dielectric saturation") as in ([34]); more precisely, in aqueous solutions, 
the compact Stern layer (defined as the inner part of the double layer not described by PB theory) 
is inferred to have an eff'ective permittivity ss smaller than that of bulk water by roughly an 
order of magnitude, e.g. reduced from = TS^o to ss = 6eo 11551 . This conclusion comes 
mostly from fitting compact layer models with nonlinear dielectric response 112271 12281 II 88L but 
also from some models assuming similar nonlinear dielectric properties in the diff'use layer. 

Grahame 12561 was perhaps the first to analyze the structure of the diff'use layer with a field- 
dependent permittivity, using PB theory and the empirical form 

where m is an empirical exponent, conjectured to be in the range < m < 2 to avoid overly sud- 
den onset of dielectric saturation above the characteristic field strength Es. For any permittivity 



model such as ( 35 ), where dielectric saturation sets in above a characteristic field strength Es, the 
importance of variable permittivity in regions of diff'use charge is measured by the dimensionless 
parameter 



(36) 

(kT/zeAo) ^ \ 2c()kT ^ c^kT 

which compares the critical electric field Es to the characteristic diff'use-layer field (at low volt- 
age), Eq = kT/zeAo, where the screening length Ad is defined using the bulk permittivity st. 
Equivalently, the dimensionless parameter E^ compares the critical electrostatic pressure for di- 
electric saturation, ps = \siyEl, to the bulk osmotic pressure, c^kT . For £^ » 1, the diff'use layer 
maintains a constant permittivity for typical voltages of order kT/ze. 

To study eff'ects of nonlinear dielectric response in the present context, we extend Grahame' s 
analysis by deriving the diff'erential capacitance of the double layer for any MPB theory with 
an arbtirary field-dependent permittivity s(E). Given the (non-Boltzmann) equilibrium charge 
density profile p(i//), we use Poisson's equation (in the normal coordinate), - (si//'y = p(if/), to 
write p{il/)il/^ = h(E)E\ where E = -if/^ and h{E) = s\E)E^ -h s(E)E. Next, we integrate from 
the bulk (i// = 0, E = 0) to the inner edge of the diff'use layer (i// = ^z), E = Ed) to obtain 

H(Ed) ^ ^ ' h(E)dE = i (s(Ed)EI + ^ " s\E)E^dE^ = p(<A)J<A - PeCi^o) (37) 

The electrostatic pressure at the surface pd^o) must be non-negative, as is the function H(Ed) 
from its definition since s > 0. Therefore, there are an even number of nonzero solutions to 
the algebraic equation H(Ed) = PeC^o), typically two for a monotonic s(E), and we select the 
physical solution by requiring that Ed and ^d have the same sign. (It can be shown that unique 
solution exists if h(E)/E = £-\- s'E > for all E.) To derive the diff'erential capacitance Cd(^d), 
we integrate Poisson's equation across the diff'use layer (i.e. apply Gauss' law) to obtain the total 
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Figure 9: The effect of dielectric saturation on the differential capacitance of the diffuse layer, given by |38j-|40j, 
for (a) PB theory (y = 0) and (b) Bikerman's MPB theory with y = 10""^. The field-dependent permittivity siE) 
given by Grahame's model {isj with m = ^ . The characteristic field strength for dielectric saturation Es is scaled to 
£^0 = kT/zeAo, where Ad = ystkljlize)^ is the screening length based on the bulk permittivity. 



charge density, q = -£{Ed)Eo- Inserting EdC^d) from (37 ), we finally obtain a general formula 
for the diff'erential capacitance of the diff'use layer, 



(38) 



valid for any p(i//) and s(E), which reduces to ( [T9]) in the case of a constant permittivity, s = Sh. 



It is instructive to consider Grahame's model ([35| with m 
solution for Co- In this case, we obtain the function 



H(E) 1 (E 



^ , since it permits an analytical 



(39) 




which can be inverted analytically by solving a quadratic equation for E^, taking a square root. 



and choosing the physical solution as described above. Substitution into (38) then yields an 
explicit formula for Cd(^d) in terms of the charge density pQi^o) and electrostatic pressure 
Pei^o) at the surface, valid for any MPB model with field-dependent permittivity £(E). For 
example, in Bikerman's model, these functions are given by 



P(^) 



2zeco smh(zei// / kT) 
l-\-2ysmh^(zei/r/2kT) 



and pe(i//) 



C()kT 



In 



l + 2vsinh^f^\ 
\2kT) 



(40) 



where PB theory corresponds to the limit v ^ 0. 

In Fig. |9] we show the eff'ects of dielectric saturation on the diff'use-layer diff'erential capac- 
itance, using our analytical formula for Bikerman's MPB theory. In all cases, the capacitance 
shows a transition from low-field dependence with s ^ to a. very similar high-field depen- 
dence with s ^ ss, which diff'ers only by a multiplicative factor ss Isb- The transition occurs 
when the electric field strength at the surface reaches Eg. Applying PB theory for low voltages, 
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this happens at a transition voltage of roughly = ze^o/kT = 2 sinh"^(£^/2), consistent with 
the plots in Fig.[9j For PB theory, shown in (a), upon increasing the voltage, dielectric saturation 
leads to a fairly sudden decrease in Co, followed by another exponential increase, corresponding 
to the PB formula ([TO]) with replaced by ss . This physical interpretation of dielectric satura- 
tion followed by "electro-constriction" is also contained in classical models of the Stern layer, 
which can lead to excellent fits of capacitance data for interfaces with little surface adsorption, 
such as NaF at a mercury drop, although only for positive polarization I257[|228lll88l . Perhaps 
the asymmetry with negative polarization not captured by PB theory with dielectric saturation, 
which is also seen in NaF on silver surfaces 112221 12231 , and may be better captured by asymmet- 
ric MPB models with diff'erent ion sizes. Here, we do not attempt such fitting, but simply discuss 
our exact result for dielectric saturation in Bikerman's MPB theory with one ion size. 

As shown in Fig. |9fb), the shape of the theoretical capacitance curves depends on whether 
dielectric saturation occurs at lower or higher voltage than ion crowding. If dielectric saturation 
sets in first, then the transition ss leads to second small peak at low voltage, separate from 
the main peak corresponding to ion crowding; if crowding sets in first, then dielectric saturation 
only causes the capacitance to drop more quickly at higher voltage. Overall, the eff'ect of dielec- 
tric saturation is similar to that of increasing the ion size. As noted below and in Ref. | 67 1, this 
eff'ect may contribute to the artificially large ions sizes inferred from Bikerman's model without 
dielectric saturation, and allow more realistic, smaller ion sizes to be used. 

It is interesting that so much structure in the diff'erential capacitance can be predicted by sim- 
ple mean-field models of the diff'use layer, without introducing a compact layer with additional 
parameters to fit large voltage behavior. Indeed, various modified models of the diff'use layer al- 
ready mimic the formation of a compact layer with increasing voltage, as ions become crowded 
in an inner condensed layer where the local permittivity drops. The "Stem plane" or "outer 
Helmholtz plane" eff'ectively expands and contracts in response to the applied voltage, in ways 
that may be captured by the same continuum equations describing the diff'use layer. Below in 
section |4j we will argue that the viscosity increases in the condensed region, and thus the "shear 
plane" also eff'ectively moves in response to the voltage. 



3. 2. Implications for nonlinear electrokinetics 

The cutoff' and eventual decrease of diff'use-layer capacitance at large voltages for blocking 
surfaces (without Faradaic reactions or adsorption of ions) is robust to variations in the model 
and has important consequences for nonlinear electrokinetics. Here, we provide two examples 
of induced-charge electrokinetic phenomena, where any MPB theory with volume constraints 
is able to correct obvious failures of PB theory. As shown above, dielectric saturation only 
enhances these eff'ects, but as a first approximation we assume constant permittivity s = st 
in our calculations. This is also consistent with theoretical estimates (63] [64| and atomistic 
simulations (1051 (at low voltages), which suggest that dielectric saturation in the diff'use layer is 
weak compared to other eff'ects. Our results suggest that incorporating crowding eff'ects into the 
electrokinetic equations may be essential in many other situations in electrolytes or ionic liquids, 
whenever the voltage or salt concentration is large. 



3.2.1. High-frequency flow reversal of AC electro-osmosis 

Steric eff'ects on the double-layer capacitance alone suffice to predict high-frequency ffow 
reversal of ACEO pumps, without invoking Faradaic reactions. Representative results are shown 
in Fig.[TO| and the reader is referred to Ref. [>67J for a more detailed study. Numerical simulations 
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Figure 10: (a) One period of an asymmetric array of planar microelectrodes in an ACEO pump studied in experi- 
ments (34][38]|40]|4T1 and simulations with the low-voltage model 1 34|[35]|42]|44 1 with Wl = 4.2 idm,W2 = 25.7 //m, 
Gl = 4.5 jum , and G2 = 15.6 jum. (b) The dimensionless flow rate versus frequency for difl'erent models. In the 
low- voltage limit V «c kT/e = 25 mV, low- voltage models predict a single peak (black dash-dot line). For a typical 
experimental voltage, V = lOOkT/e = 2.5 V, PB theory breaks downs and its capacitance flO) shifts the flow to low 
frequency (red dashed line) and Stern capacitance is needed to prevent the capacitance from diverging. Accounting for 
steric efl'ects J25j with y = 0.01 (solid blue line) reduces the shift and predicts high frequency flow reversal, similar to 
experiments 138 41 1 . 



of a well studied planar pump geometry 1341 [381 |401|4Tl with the Standard Model in the linearized 
low- voltage regime |11, 12,|35l|44]| predicts a single peak in flow rate versus frequency at all 
voltages. If the nonlinear PB capacitance ([TO]) is included 1421 |43l, then the peak is reduced and 
shifts to much lower frequency (contrary to experiments), due to slower charging dynamics at 
large voltage [ 129111311 . As shown in Fig. 10 the BF capacitance for Bikerman's MPB model 
of steric efl'ects ( [25] ) reduces the peak shift and introduces flow reversal similar to experiments. 

This result is the first, and to our knowledge the only, theoretical prediction of high-frequency 
flow reversal in ACEO. The physical mechanism for flow reversal in our model can be easily un- 
derstood as follows: At low voltages, the pumping direction is set by the larger electrode, which 
overcomes a weaker reverse flow driven by the smaller electrode. At large voltages, however, the 
more highly charged, smaller electrode has its RC charging time reduced by steric efl'ects, so at 
high frequency it is able to charge more fully in a single AC period and thus pump harder than 
the larger electrode. 

As shown in Fig.[TT] the MPB model is able to reproduce experimental data for ACEO pump- 
ing of dilute KCl rather well, including the dependence on both voltage and frequency. Through 



Fig. 11 we compare simulations to experimental data at two difl'erent concentrations. In the 
left column we show experimental data and on the right we show the corresponding simulations 
using Bikerman's MPB theory for the double later capacitance. As in experiments, the flow re- 
versal arises at 10-100 kHz frequency and high voltage, without shifting appreciably the main 
peak below 10 kHz frequency (which is hard to see in experiments at high voltage due to elec- 
trolysis). This is all the more remarkable, since the model has only one fitting parameter, the 



eff'ective ion size a, and does not include any additional Stern-layer capacitance. As seen in 1 1 
(a) and (b), the magnitude of the flow is over-estimated by a roughly a factor of two (A ^ 0.4), 
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Figure 11: (a) Velocity of ACEO pumping in 0.1 mM KCl by a planar electrode array around a micro fluidic loop versus 
frequency at different voltages from the experiments of Studer et al. | 38 1. (b) Simulations by Storey et al. |^ of the 
same flow using the Standard Model with Bikerman's MPB theory \25) for the double-layer difl'erential capacitance with 
only one fitting parameter, a = 4.4 nm or y = 0.01. Countour plots of ACEO pumping velocity contours in frequency- 
voltage space for (c) 0.1 mM and (e) 1.0 mM KCl from experiments of Studer et al. 1 38 1, compared to simulations under 
the same conditions using Bikerman's MPB theory with v = 0.01 in (d) and (f), respectively. Red indicates forward flow 
and blue reverse flow. The solid contour lines show positive velocity contour and the dashed show reverse flow. The 
heavy solid contour is the zero velocity contour in the simulations. 
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but this is much better than in most predictions of the Standard Model (in Table 1), which fail to 
predict flow reversal under any circumstances. 

In spite of this success, we are do not claim a complete understanding of flow reversal in 
ACEO. One difficulty with these results is that the effective ion size in Bikerman's model needed 
to fit the data is unrealistically large. For the simulations to reproduce the experiments we seem 
to typically require v = 0.001 - 0.01, which implies an overly small bulk ion spacing /q = 
(2co)~^^^ = v'^^^a, or overly large ion size a. For example, for the cq = 0.1 mM KCl data shown 
in Fig. 1 1 (a) and (b), we use a = 4.4 nm in the model, which is clearly unphysical. As noted 
above and shown in Fig.|7j this can be attributed at least in part to the signficant under-estimation 
of steric eff'ects in a liquid by the simple lattice approximation behind Bikerman's model. 

Indeed, hard- sphere liquid models tend to improve the agreement betweeen simulation and ex- 
periment, and this increases our confidence in the physical mechanism of ion crowding at large 
voltage. Using the CS MPB model for monodisperse charged hard spheres in the same simula- 
tions of ACEO pumping allows a smaller value of the ion size. For example, the 0.1 mM KCl 



shown in Fig. 1 1 (a) can be fit by using a = 2.2 nm (instead of 4.4 nm for Bikerman), and 



the magnitude of the velocity also gets closer to the experimental data (A ^ 0.7). Assuming a 
reduced permittivity in the condensed layer could further yield a ^ 1 nm (^10 atomic diam- 
eters) |67|. This value is more realistic but still considerably larger than the bulk hydrated ion 
sizes in KCl and NaCl. For the commonly used electrolytes in ICEO experiments (see Table [T]), 
the cation-anion radial distribution function from neutron scattering exhibits a sharp hard-sphere- 
like first peak at 3 A, although the water structure is strongly perturbed out to the second neighbor 
shell (up to Inm), as if under electrostriction. Anion-anion correlations are longer ranged and 
softer, with peaks at 5 Aand 7 A, but unfortunately such data is not available for crowded like 
charges within the double layer at high voltage. Perhaps under such conditions the effective 
hard-sphere radius grows due to strong correlations, beyond the mean-field approximation. 

We have already noted that similar over- sized ionic sizes are needed to fit electrochemical ca- 
pacitance data with MPB models. (See Fig. [8j where a = l.l nm in Bikerman's model yields 
an excellent fit of capacitance data for KPFg | Ag.) In the case of hard- sphere MPB models, 
DiCaprio et al. have likewise concluded that, in spite of good qualitative predictions, effective 
sphere radii over 1.2 nm are required to fit capacitance data 11921 . The fact we reach similar con- 
clusions in the completely different context of ACEO pumping suggests that various neglected 
effects, in addition to dielectric saturation (section 3.1.6 ), may be extending the apparent scale 
for crowding effects in MPB models, as discussed below. 

Along with overly large ion sizes, another difficulty is that the simulated results in Fig.[TT](d) 
and (f), are reproducing data from two different concentrations using the same value of v = 0.01 
in Bikerman's MPB theory. Since the concentration varies between the two experiments by a 
factor of 10, the simulations are using different values of the ion size, a: 4.4 nm at 0.1 mM and 2 
nm at 1 mM. Physically, it would make sense to use the same ion size in the model regardless of 
concentration. If we use the same ion size to fit data from both concentrations, the comparison 
worsens significantly, representing one of the key difficulties of fitting model to experimental 
data. 

The simple steric models (without dielectric saturation) can thus provide good qualitative 
agreement between the measured and predicted frequency responses. However, they can not 
predict the proper frequency response as we change concentration and hold the ion size fixed. In 
the model calculations, the frequency response at high voltage always shows regions of forward 
and reverse flow. As concentration changes, features in the modeled frequency response (e.g. 
forward and reverse peaks, crossover frequency) shift along the frequency axis as the RC charg- 
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ing time changes with concentration. Since the capacitance of the double layer is insensitive to 
bulk concentration at high voltage (see Fig. [6]), the RC charging time depends on concentration 
through the bulk resistance. Thus the location of peaks on the frequency response should vary 
approximately linearly with concentration; experiments at high voltage indicate a much weaker 
dependence. The difficulty in matching the frequency response of experiments to simulations 
may be due to the neglect of Faradaic reactions and other effects discussed below. 

At least at a qualitative level, changes in the double-layer charging time due to crowding 
effects likely also play a role in the sensitivity of ACEO flow to the solution chemistry. For 
example, flow reversal in our models is related to an ion-specific scale for crowding effects. In 
cases of asymmetric electrolytes, there may also be two different frequency responses at large 
voltages, one for positive and another for negative charging of the double layer. Perhaps this 
effect is responsible for the double-peaked frequency spectrum of ACEO pumping in water with 
non-planar electrodes at high voltage In multicomponent electrolytes, the situation is even 
more complicated, since large voltages can induce segregation of different counterions, oppo- 
site to PB predictions, e.g. with smaller ions condensing closest to the surface, even if the 
larger ions have high bulk concentration or carry more charge. These effects can be predicted 
by multi-component hard-sphere MPB models 1621 11921 11931 consistent with x-ray reflectivity 
measurements on mixed double layers 1258112591 , so there is hope that applying such models in 
our weakly nonlinear formalism for ACEO may also be fruitful. 

3.2.2. Field- dependent electrophoretic mobility 

In the classical theory of electrophoresis |[8l[6l, the electrophoretic mobility bep of a homoge- 
neous particle with thin double layers is a material constant, given by Smoluchowski's formula, 

bep=b=^-^. (41) 

m 

In particular, the electrophoretic mobility does not depend on the background field Eij or the 
shape or size of the particle. These are partly consequences of the assumption of fixed surface 
charge, or constant zeta potential. For polarizable particles, the theory must be modified to 
account for ICEO flows |19|, which produce a shape- sensitive ICEP velocity scaling as oc 
£hREl/r]iy, where R is the particle size |[T5l[23l|24l. Transverse ICEP of metallo-dielectric Janus 
particles in AC fields has recently been observed in experiments |25 1 up to fairly large induced 
double-layer voltages E^R ^ l5kT/e. 

ICEP theories aimed at AC fields tend to assume zero total charge, but ICEO flows can also 
alter the DC electrophoretic mobility of a charged, polarizable particle. In the limit of weak fields 
Eh <$c kT/eR, A. S. Dukhin first showed that an ideally polarizable sphere with equilibrium zeta 
potential and radius R has a field-dependent electrophoretic mobility, 

T]b \ 8 Cd(^o) 



bep(Et) ~ - ko - -^^,(EtR)' + . . . I (42) 



if the diffuse-layer differential capacitance is voltage dependent |65|. (Note that field-depedent 
mobility is a general phenomenon |66 | that can also arise for fixed-charge particles due to surface 
conduction |111| or convection |110).) This general correction has only been applied in the 
context of PB theory, Eq. ([TO]), which predicts decreased mobility, Abep < since dCo/di// > 
for ^0 > 0- It has also recently been derived as the small field limit of a general PB analysis for 
thin double layers by Yariv LI 45 J . 
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Figure 12: Field-depenent electrophoretic velocity U of an ideally polarizable, charged sphere of radius R with thin 
double layers in a background field E},. (a) In small fields, the mobility bep = U/Et> is set by the uniformly distributed 
double-layer charge, (b) In large fields, Eh » kT/eR, the dipolar induced charge overwhelms the pre-existing charge and 
alters bep, if cations and anions do not condense at the same density and must redistribute to conserve total charge, (c) 
In PB theory, the unphysical collapse of point-like ions to the surface causes exponential decay of bep(Eb) via Eq. |44|; 
finite-size eff'ects in Bikerman's MPB model (Fig. |5} prevent this decay and lead to the opposite trend: increase of 
mobility in large fields via Eq. J47}. 
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The basic physics of this nonlinear effect is illustrated in Fig. 12 a-b). If the double-layer 
voltage varies enough to cause spatial variations in its differential capacitance, then counterions 
aggregate with varying density (per area) around the surface of the particle upon polarization 
by the applied field, and this nonlinearity breaks symmetry in polarity with respect to the mean 
voltage. For example, if the positively charged part of the diffuse layer (relative to the mean 
charge) is less dense (e.g. due to larger or less charged cations than anions), it will cover more 
of the surface than the negatively charged part; cations are then more likely to dominate in 
regions of large tangential field near the equator and thus make an enhanced contribution to the 
electrophoretic mobility of the particle, regardless of its true surface charge. Other effects can 
also be important (see below), but this one is particularly sensitive to MPB models for the double 
layer. 



Dukhin's formula (42) can be derived from the general weakly nonlinear formalism of Refs. 
|[T9l [2TI for ideally polarizable particles with thin double layers (yielding the same result as 
Ref. (65 1). In the low-frequency or DC limit, the background field £^ causes nonuniform polar- 
ization of the double layer around the particle to screen the bulk electric field E = -V0, which 
thus solves Laplace's equation = with the effective boundary condition, h - E = 0. If we 
let 00 denote the induced potential of the particle, relative to the background applied potential, 
then A0(r) = (pQ - (p(r) is the non-uniform voltage across the double layer, which enters the 
electro-osmotic slip formula, either the HS formula ^ or one of its generalizations below. For a 
sphere with HS slip, the electrophoretic mobility is simply bep = st(po/T]t. 

The crucial step is to determine the particle's potential (po, after polarization of the double 
layers, which generally differs from its initial equilibrium value, ^0, due to nonlinearity of the 
charging process. (This phenomenon, first noted by Dukhin [65 1, was overlooked in recent ICEO 
papers (TSl [161 EH EH-) Using the mathematical formalism of Ref. 1231 , the potential of the 
particle must adjust to maintain the same total charge Q, which can be related to the differential 
capacitance of the double layer as follows: 

G = j)N CD(i/f)di/f\dA = CD(i/f)dHdA, (43) 

for a given initial zeta potential fo, bulk polarization 0(r), and (total) double layer differential 
capacitance, Cd(i//)- For simplicity, we assume the diffuse layer carries all the double-layer 
voltage, but a compact Stern layer can be easily included in Eq. ( [43] ) by replacing A0 with 
A0/(1 -h S). Regardless of the particle shape, the assumption of a uniform, constant Co implies 
(po = ^0 = Q/(CdA) (where A is the surface area), and thus no impact of polarization on the 
electrophoretic mobility in the case of a sphere |[23ll . With a nonlinear differential capacitance. 



however. Equation (43 ) is a nonlinear algebraic equation for and thus bep, in terms of and 



El,. In the geometry of a sphere, Dukhin's formula (42 ) can be derived by asymptotic analysis in 
the limit of small fields, «: kT/eR for any choice of CdW IIH, and for some models exact 
solutions and the large-field limit can also be derived |260 |. 

Using this mathematical formalism with our MPB models, we predict that steric effects in the 
electrolyte can signficantly influence the mobility of polarizable particles in large applied fields 
and/or highly concentrated solutions. Here, we focus on new qualitative phenomena predicted by 
the theory. (More details can be found in Ref. |260|.) From Fig. [5] and Eq. ( [25] ), we see that the 
mobility of a highly charged particle |^ol ^ kT/e can increase with the field squared in Dukhin's 
formula ( [42] ) since dCo/dip < 0, which is the opposite prediction of PB theory. The mobility is 
also clearly sensitive to the ionic species through Co- 
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The discrepancy with PB theory becomes more dramatic in a large appHed field, Eiy » kT/eR, 
even if the particle is not highly charged |^ol ~ kT/e. Previous authors have only considered weak 
fields (661 EH, so it was apparently not noticed until Refs. 111451 12601 that PB theory ([TO]) leads 



to a surprising prediction, shown in Fig. 12 c): The mobility of a charged, ideally polarizable. 



spherical particle vanishes exponentially in the limit » kT/eR (=100 V/cm for R = 2.5yum), 

Z,™ ~ ^ sinh (^)EtR .-^-^-.^/^^ (44) 

which is the large-voltage limit of an exact solution for the ICEP mobility of an ideally polariz- 
able sphere with thin double layers in PB theory. 



ze 



3zeEtR sinh(z^f / 2^r) 



4kT smh(3zeEtR/4kT) 



(45) 



The mechanism for this seemingly unphysical eff'ect is the massive overcharging of the diff'use 
layer in PB theory at large voltages in Fig.[5jb), which causes the anti- symmetric induced charge 
(not causing motion) to overwhelm the symmetric pre-existing charge (giving rise to mobility). 
We view this prediction as another failure of PB theory, since we are not aware of any evidence 
that polarizable particles lose their electrophoretic mobility so strongly in such fields, which 
are routinely applied in electrophoresis experiments. (Note that kT/eR = 100 V/cm for a 5yum 
diameter particle at room temperature.) The PB prediction of vanishing mobility is closely tied 
to the unphysical pile-up of point-like ions near a highly charged surface in PB theory. 

Indeed, the situation is completely diff'erent and more physically reasonable in any mean-field 
theory with finite-sized ions, regardless of the model. Using our general MPB formula ([19]), the 
mobility in large fields E^ » kT/eR can be expressed as 

3sbEbRq(^o) 

bep -j- r (46) 

2r]bq[^EtR) 



and typically grows as oc yS^. For example, in Bikerman's model the BF formula ( 25 ) yields 



SvkTEbR 

As shown in Fig. [T2tc), the mobility b^p rapidly decreases in weak electric fields until PB theory 
breaks down, and men gradually increases in larger fields | 260 |. Physically, steric eff'ects prevent 
overcharging of the double layer in PB theory, thus preserving the asymmetry of the initial charge 
distribution. 

A related general consequence is that asymmetry in the electrolyte, e.g. ions of diff'erent 
eff'ective sizes, can aff'ect a particle's electrophoretic mobility. Remarkably, an uncharged particle 
can acquire a nonzero mobility in an asymmetric electrolyte in a large applied field (fo = but 
bep ^ 0). An asymmetric field-dependence of the double-layer capacitance is enough to predict 
such exotic eff'ects on very general grounds. 

As a first approximation for the weakly nonlinear regime using MPB theory, we simply postu- 
late two diff'erent (homogeneous) ion sizes a+ for positive and negative polarization of the diff'use 
layer, 

CdCVd) - I ^ ^ < (48) 
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Figure 13: Electrophoretic mobility bep of an ideally polarizable sphere of total charge Q in an asymmetric binary z : z 
electrolyte, scaled to b^p = EbkT/zerjb, in the weakly nonlinear limit of thin double layers. The case of an uncharged 
particle Q = Ois compared to those of total charge Q = ±Q^ where Q'^ = SbkT/zeAo. Using the approximation jJs}, the 
mobility is calculated with Bikerman's MPB model |25) for an effective volume ratio (a-/a+)^ = 10 in two cases: (a) 
high salt concentration with y_ = 0.1 and y+ = 0.01 and (b) moderate salt concentration with y_ = 10~^ and v+ = 10~^. 
At small electric fields and/or low salt concentrations, the size asymmetry is irrelevant, and the predictions of PB theory 
from Fig[T2]are apparent; at large fields and/or low concentrations, the particle acquires an apparent positive charge, due 
to the covering of more of the particle's area by the larger cations, regardless of its true surface charge; 



where y- = 2a+co (or O- for hard sphere models). This may seem hke a crude approximation, but 
it is quite accurate for a binary electrolyte at large voltages since the diffuse layer mostly contains 
counterions of one sign; at low voltages, where all species are present in a dilute mixture, the 



ion sizes play no role. Inserting ( [48) ) into ( [43] ) yields the results shown in Figure [13] for the 
(weakly nonlinear) mobility of a charged, ideally polarizable sphere in a DC field. Interestingly, 
at large voltages, positive or negative, the mobility tends to a \E\ scaling, set by the ratio a+/a-, 
independent of the total charge of the particle Q. We also see that an uncharged particle with 
2 = can still have a nonzero mobility, once nonlinear charging of the double layers sets in. 

We stress that all the calculations of this paper consider the weakly nonlinear limit of thin dou- 
ble layers, where the bulk concentration remains uniform and surface conduction is neglected. 
As such, the trends we predict for PB and MPB models are meaningful at moderate voltages, but 
may need to be significantly modified at large voltages for strongly nonlinear dynamics. More- 
over, even in the weakly nonlinear regime, we will now argue that at least one more physical 
important eff'ect may need to be considered. 



4. Viscoelectric effect in a concentrated solution 

4.1. Mean-field local-density theories 

4.1.1. Modified Helmholtz-Smoluchowski slip formulae 

There is a considerable literature on electrokinetic phenomena at highly charged surfaces with 
large (but constant) zeta potential 161 [5l 11121 II 131 [TOll . In this context, it is well known that the 
linear electrophoretic mobility departs from Smoluchowski's formula at large surface potentials, 
^ > and decreases at large voltage, due to eff'ects of surface conduction (large Dukhin num- 
ber). Using PB theory for thin double layers, Dukhin and Semenikhin LI 12 1 derived a formula 
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for the electrophoretic mobility of a highly charged, non-polarizable sphere, which was famously 
verified by O'Brien and White | 261 1 via numerical solutions of the full electrokinetic equations 
for a dilute solution. This established the mathematical validity of the formula, but in this article 
we are questioning the physical validity of the underlying equations at large induced voltages 
and/or high salt concentrations. 

As described in section |2j recent experiments on induced-charge electrokinetic phenomena 
reveal a strong decay of electro-osmotic mobility with increasing salt concentration at highly 
charged surfaces, which cannot be explained by the standard model, based on the HS formula 
([5]), even if corrected for strongly nonlinear eff'ects perturbing the salt concentration. Continuum 
electrohydrodynamics, however, does not require the HS formula, even for thin double layers, 
but instead provides a general expression for the electro-osmotic mobility [6], 



as an integral over the potential diff'erence ^ entering the diff'use layer from the bulk. This allows 
us to derive "modified Helmholtz-Smoluchowski" (MHS) formulae for b(^D, cq, . . .) based on 
general microscopic electrokinetic equations. 

The standard way to interpret electrokinetic measurements is in terms of the eff'ective zeta 
potential. 



but we view this as simply a measure of flow generated in units of voltage, and not a physically 
meaningful electrostatic potential. One can also use PB theory for the hypothetical "mobile" part 
of the double layer to express the mobility in terms of an efl'ective "electrokinetic charge" |[T0l[6l, 
e.g. using the Gouy-Chapman solution for a symmetric binary electrolyte. 



which measures the observed flow in units of charge. 

It is well known that the classical theory tends to overpredict experimentally inferred zeta 
potentials (\^eff\ < \^d\ or \qek\ < \q\), and the discrepancy is often interpretted in terms of a 
"slip plane" or "shear plane" separated from the surface by a molecular distance d > 0, where 
the no-slip boundary condition is applied 12621 [TOl . If the low- voltage theory with s = and 
T] = T]iy is applied beyond this point, then ^eff acquires physical meaning, as the potential of the 
slip plane relative to the bulk. Behind the slip plane, the fluid is assumed to be "stagnant" with 
efl'ectively infinite viscosity, although it may still have finite ionic surface conductivity 1101 . 

As in the case of eff'ective hydrodynamic slip (d < 0) over hydrophobic surfaces 12631 12641 
the slip-plane concept, although useful, obscures the true physics of the interface. In particular, it 
has limited applicability to nonlinear electrokinetic phenomena, where the double layer responds 
non-uniformly to a large, time-dependent applied voltage. We have already argued this leads to 
signficant changes in the structure of the double layer, and clearly it must also have an impact on 
its rheology. Without microscopic models of how local physical properties, such as viscosity and 
permittivity, depend on local variables, such as ion densities or electric field, it is impossible to 
predict how the eff'ective electro-osmotic slip depends on the global double-layer voltage or bulk 
salt concentration. This is also true in complicated geometries, such as nanochannels or porous 
structures, where the concept of a flat "slip plane" is not realistic. 




(49) 




(51) 
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4.1.2. The viscoelectric effect 



Based on Eq. (49 ), there are good reasons to expect reduced mobility at large voltage, |^z)| > 



compared to the HS formula, ^^ff = based only on field-dependent properties of a polar 



solvent |6|. As discussed in section 3.1.6 large normal electric fields in a highly charged double 
layer can decrease the permittivity by aligning the solvent dipoles ("dielectric saturation"). The 
local viscosity can also increase 12651 1251 L through viscoelectric thickening of a dipolar liquid 
with polarization transverse to the shear direction. 

Long ago, Lyklema and Overbeek L63| ,64J (LO) proposed the first (and to date, perhaps the 
only) microscopic electro-rheological model leading to an MHS slip formula. They focused on 
the viscoelectric eff'ect in water, which they estimated to be more important than dielectric satura- 
tion for electrokinetics (consistent with recent atomistic simulations 1 105 1, albeit at low voltage). 
Physically, they argued that, for typical electric fields in aqueous double layers, dielectric sat- 
uration interferes more with the strongly cooperative rearrangements involved in viscous flow 
than the weakly correlated dipole alignments contributing to reduced permittivity. To model the 
viscoelectric efl'ect, LO assumed a field- squared viscosity increase in PB theory, 

77 = 77,(1+/^^), (52) 



and were able to integrate ( 49 ) to obtain an analytical (although cumbersome) MHS formula. 
Physically, their model predicts that b saturates to a constant value at large which decays 
with increasing cq. The saturation in the LO model, however, is tied to the unphysical divergence 
of the counterion concentration (and thus E) in PB theory, and thus should be revisited with 
volume constraints and other eff'ects. 

It is straightforward to use the LO postulate ([52]) for the viscoelectric eff'ect in our modified PB 



models to obtain corresponding MHS slip relations, although the integration of (49) cannot be 



done analytically. Numerical results for a typical case are shown in Fig. [14] using the value / = 
^q-15j^2y-2 suggested by LO for water |[63j[64l. It is interesting to note that this choice makes 
all three models of double layer structure, PB, Bikerman MPB, and CS MPB, yield very similiar 
electro-osmotic mobility versus voltage, in spite of completely diff'erent ion density profiles. The 
reason is that the viscoelectric eff'ect sets in so quickly with increasing voltage that the shear plane 
is eff'ectively still in the dilute, outer part of the diff'use layer, where all theories reduce to PB. 
Indeed, as shown in the figure, more diff'erences become apparent we if choose a smaller value, 
/ = 10"^^m^V"^. The reduction in normal electric field due to crowding-induced expansion of 
the inner diff'use layer leads to slower saturation of ^eff compared to PB theory, especially in the 
CS MPB model, since it has stronger steric eff'ects than the Bikerman MPB model. 

What this exercise shows is that the inner part of a highly charged double layer is eff'ectively 
frozen by field-dependent viscosity in a similar way regardless of the model for the diff'use- 
charge profile. The LO model attributes this eff'ect entirely to the solvent (water), independent 
of the local diff'use charge density or ionic currents, but this hypothesis does not seem entirely 
satisfactory. An implication is that water is eff'ectively immobilized within a few molecular 
layers of the surface near a charged surface, even if there is no added salt. This eff'ect seems 
to contradict recent experimental and theoretical literature on hydrodynamic slip ||263[ [TJ 12661 
12641 which has shown that shear flow in water persists down to the atomic scale with a no- slip 
boundary condition for smooth hydrophilic surfaces and with significant slip lengths (up to tens 
of nanometers) on hydrophobic surfaces 12671 12681 12691 . If the viscosity is a property of pure 
water (as opposed to the local electrolytic solution), then the field dependence should also be 
observable in the bulk, but we are not aware of any experiments or simulations showing that 
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Figure 14: Modified Helmholtz-Smoluchowski (MHS) slip formulae assuming a charge-independent viscoelectric effect 
in the polar solvent {52) . This example assumes a bulk concentration cq = 1 mM of z : z electrolyte using different 
models of double-layer structure. The Lyklema-Overbeek (LO) model, based on Poisson-Boltzmann (PB) theory of 
point-like ions, is compared to MHS slip with the Bikerman (MPB) and Carnahan-Starling (CS) modified PB theories. 



using an effective ion size of « = 4 A. The viscoelectric coefficient is set to the value / = 10 
for water as well as a smaller value / = 10~^^m^V~^, which reduces the viscoelectric effect. 



V suggested by LO 



bulk de-ionized water becomes rigid in fields larger than /"^ = 30V/yum. Of course, it is hard to 
apply such fields in the bulk at low frequency, due to capacitive screening by the double layers 
and Faradaic reactions, but the LO model is time-independent and should also apply at high 
frequency where these eff'ects are reduced. 

Molecular dynamics (MD) simulations of electrokinetic phenomena in nanochannels gen- 
erally imply a local viscosity increase close to a charged surface, although its origin and 
mathematical description are not well understood. An early MD simulation of Lyklema et 



al. 110411 showed a stagnant monolayer of water, which could pass ions freely, with a sur- 
face conductivity comparable to the bulk, but subsequent MD simulations of electro-osmosis 
have shown motion of the liquid (both ions and solvent) down to the atomic scale near the 
wall 1^ ^ El ElQl ^ una, albelt with an apparent viscosity often smoothly 
increasing across the closest molecular layers. We are not aware of any MD simulations of 
electro-osmosis in a large transverse voltage (» kT je), but it seems that the presence of crowded 
counterions near the surface (described in section [3]) must aff'ect the apparent viscosity, since the 
liquid no longer resembles the pure bulk solvent. 

4.1.3. Charge-induced thickening 



A natural physical hypothesis, sketched in Figure [15] is that the crowding of counterions in a 
highly charged diff'use layer increases the local viscosity, not only through the viscoelectric eff'ect 
of the bare solvent, but also through the presence of a large volume fraction of like-charged ions 



compressed against the surface. For now, we neglect explicit field dependence, such as (52), 
to focus on eflfects of large charge density. In a very crude attempt at a model, we consider an 
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electrolyte with solvated ions of finite size and postulate that sir] diverges as a power law, 



(53) 



as the local charge density approaches a critical value, p*, which generally must depend on the 
polarity sign(p) = ± in an asymmetric electrolyte. 

A natural choice is to postulate diverging viscosity ("jamming") at the steric limit, pj = p*^^ = 
k±\^^max- that case, similar exponents a and jS controlling the singularity also arise in the 
rheology of dense granular materials II271L If we also assume a = /3 = 1, then there are no new 
fitting parameters, 

^ = ^(1-^) (54) 

since the steric constraints p^^^ come from the MPB model. In our original letter 1571 . we 
considered only this postulate with Bikerman's model for steric eff'ects and thus assumed an 
even simpler form. 

The resulting electrokinetic model is extremely simple in that it only involves one parameter 
beyond dilute- solution theory, the eff'ective ion size, a. As such, it lacks flexibility to fit multi- 
ple sets of experimental data, but we will use it in our analytical calculations below to further 
understand the general consequences of steric constraints. 



The general model (53) also allows for other types of behavior. With p* > p^ax^ there is a 
finite maximum viscosity in the condensed layer, or eff'ectively some flow behind the slip plane. 
With p'j < p^ax^ the model postulates flow arrest at high charge density before the close-packing 
limit is reached (and t] = 00 for \p\ > pj). The new parameters pj, a and allow some flexibility 
to fit experiments or simulations, in addition to the ion sizes a+ from the MPB models above. 



The arbitrary choice (53 ) is motivated by a number of possible physical eff'ects: 



• Jamming against a surface. In a bulk colloid |7|, the maximum density corresponds to 
random close packing at the jamming point L272J . where the shear modulus becomes fi- 
nite 1 273 1, diff'usivity vanishes 12741 12751 12761 and the viscosity diverges 12711 . Molecular 
dynamics simulations of soft disks in a periodic box have recently established a viscos- 



ity divergence of the form (53) with a = I and j3 = 1.7 l>271J . In an electrolyte, strong 
electrostatic compression of solvated counterions against a (typically rough) surface may 
cause some transient local jamming of the condensed layer, and thus increased viscosity 
(and decreased mobility of individual ions) at high charge density. 

Electrostatic correlations. Condensed counterions at large voltages ressemble a Wigner 
crystal 1277112781 (or glass) of like charges. Discrete Coulomb interactions I TSTl [2791 fT89l 
may contribute to increased viscosity, e.g. through the attraction between a displaced ion 



and its "correlation hole", which eff'ectively carries an opposite charge. (See Fig. [16]) We 
are not aware of attempts to predict the rheological response of sheared Wigner crystals 
or glasses, let alone discrete, correlated counterion layers, but we expect that electrostatic 
correlations will generally contribute to charge-induced thickening in electrolytes and ionic 
liquids. 

44 




Figure 15: Sketches of finite-sized hydrated ions near a polarizable surface as in Fig.|4] showing the solution velocity u 
profile in response to a tangential electric field, (a) At small induced voltages, the no-slip boundary condition holds at the 
surface, and the effective electro-osmotic slip Us builds up exponentially across the diffuse layer, (b) At large voltages, 
crowding of hydrated ions increases the viscosity of the condensed layer, and the apparent slip plane "SP" (dashed line) 
moves away from the surface with increasing voltage. 
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Figure 16: Sketch of a physical mechanism for charge-induced thickening due to electrostatic correlations. A condensed 
layer of solvated counterions (spheres) is confined against a highly charged surface by the normal electric field. Shear 
of the fluid (arrow) causes an ion to leave its local equilibrium position in a Wigner-like crystal of like charges, but its 
motion is inhibited by a strong attraction back to its oppositely charged "correlation hole" (dashed). 



• Solvent effects. Several molecular dynamics simulations of linear electro-osmosis (at low 
voltages) have infered 5x greater viscosity within 1 nm of a flat surface 11051 11061 , 
which grows with surface charge |270| and surface roughness |280| (but decreases with 
hydrophobicity 11081 ). Dielectric saturation near a highly charged surface (section 3.1.6| ) 
reduces the electro-osmotic mobility by lowering s, but the aligning of dipoles leads to 
collective interactions that can also increase rj, i.e. the classical viscoelectric efl'ect, at even 
lower fields 1631 . In the presence of crowded, like-charged counter-ions, frustrated solvation 
shells and confined hydrogen-bond networks in water may further increase //. The latter 
eff'ect may be related to local electrostriction, which compresses the second hydration shell 
around certain cations (K^ and Na^, also used in ICEO experiments. Table [T]) 11821 . 

Of course, this electrohydrodynamic model is rather simple. In general, we expect that s and 
7] will depend independently on both the local field E and the solution composition. We have 
already mentioned the diff'erent physical eff'ects aff'ecting s and // in a pure dipolar solvent, and 
the situation only becomes more complicated with large ion densities. In order to model general 
phenomena such as diff'usio-osmosis or surface conduction, it is necessary to provide separate 
functional forms for s and rj, since these variables no longer only appear as the ratio s/rj. In that 
case, we would propose viewing ( [53] ) as a model for rj with s = Sh fixed, since molecular dynam- 
ics simulations of electro-osmosis predict smaller changes in permittivity than in viscosity near 
a charged surface, even within a few molecular diameters 1 1 05 , 106] . This is also analytically 
convenient because variable £ is a major complication in MPB models, although progress can be 
made as shown in section 3.1.6 Finally, our neglect of explicit field-dependence such as (52), 
in favor of charge-density dependence (which is more closely related to our modeling of crowd- 
ing eff'ects), does not aff'ect the modeling of quasi-equilibrium double layers, since the dominant 
normal field can be expressed in terms of the charge-density, and vice versa, in any MPB model. 
In more general models, however, it may be important to also include a field-dependent visco- 
electric eff'ect. 

Our basic picture of charge induced thickening is consistent with prior models for the increase 
of viscosity with salt concentration in neutral bulk solutions [ 2651 12511 12531 . Traditionally, 
the viscosity of electrolytes has been described empirically using the Dole- Jones correlations 
12811 12821 1283 1 but more recent theories based on molecular models take a more fundamental 
approach. Like our hypothesis, these models also postulate infinite viscosity as the hard sphere 
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solvated volume fraction approaches unity f253l . In a neutral bulk electrolyte, such high volume 
fractions are never reached due to limits of solubility, though a factor of ten increase in viscosity 
has been observed for some systems at high salt concentration II284L The effective hard-sphere 
diameters (see Table [2]) inferred from these viscosity models may also have relevance for our 
models, although there is an important difference: We neglect crowding effects in a neutral bulk 
solution and focus instead on crowding of like-charged counter-ions. As discussed above, at 
high charge density in a large electric field, various physical effects could significantly increase 
the viscosity compared to what would be observed in a neutral solution at the same concentra- 
tion. Moreover, solubility limits are not relevant for like-charged ions, so higher concentrations 
approaching packing limits can more easily be reached, as we have postulated in response to a 
large voltage. 

Regardless of the specific form of our model, its key feature is that the concentrated solution 
of counterions near a highly charged surface is thickened with voltage, compared to the bulk, 
although it can still flow slowly. This picture is consistent with molecular dynamics simulations 
of electro-osmosis II 1051 11061 [270l[Tolll280ll285L which do not observe a truly stagnant layer, 
even at low surface potentials. As sketched in Fig. 15 from a macroscopic point of view, there is 
an apparent slip plane separated from the surface, which depends on the ion density profiles and 
thus external conditions of voltage and bulk solution composition. At low voltage, near the point 
of zero charge, there is no change in the local viscosity, aside from any pure solvent interactions 
with the surface, e.g. due to hydrophobic or hydrophilic effects (which we neglect). At high 
voltage, the crowding of counterions leads to thickening and an apparent movement of the slip 
plane away from the surface. Without a microscopic model, it would not be possible to predict 
this dependence or derive a functional form for fitting. 

The picture of a thin layer near the surface with different electro-rheological properties also 
appears in models of electrokinetics with porous or soft surfaces | 286, 287 ], which build on the 
Zukowski-Saville 1 288, 289] and Mangelsdorf- White 1290,180, 181| theories of the "dynamical 
Stern layer". In these models, ions are allowed to move within a flat Stern monolayer, while the 
diffuse layer is described by the standard electrokinetic model. Recently, Lopez-Garcia, Grosse 
and Horno have extended these models to allow for some fluid flow in a thicker dynamical sur- 
face layer, and this allows more flexibility in fitting (linear) electrical and electrokinetic measure- 
ments 12911 12921 . In the case of a flat surface with an equilibrium double layer, this is similar to 
our picture, if the surface layer is ascribed a higher viscosity, but in our model there is no need 
to postulate a sharp plane where properties of the solution change. Instead, by modifying the 
electrokinetic equations everywhere in the solution, its electro-rheology changes continuously as 
a function of local contiuum variables. 

The concept of charge-induced viscosity increase may be widely applicable in nano-scale 
modeling of electrolytes, in conjunction with modified theories of the ions densities. The general 
modified electrokinetic equations are summarized in section |5] Now, let us consider generic 
consequences of this hypothesis for electrokinetic phenomena with thin double layers and make 
some explicit calculations using Eq. ([53 ). 



4. 2. Implications for nonlinear electrokinetics 

4.2.1. Electro -osmotic slip at large voltages in concentrated solutions 

To describe electro-osmotic flows with thin double layers, we can use our microscopic models 
to derive modified Helmholtz-Smoluchowski (MHS) slip formulae, which depend nonlinearly 
on double-layer voltage (or surface charge) and bulk salt concentration. Such predictions could 
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Figure 17: MHS slip formula p9) using Bikerman's model with ion size a = 4Aand charge-induced thickening J53j 
with a=l=fi=lsit different bulk concentrations, cq = 1/jM, 1 mM, and IM. The viscosity is postulated to diverge 
either at a mean ion spacing aj = a (solid curves), in which case the condensed layer in Fig.jSja) is effectively rigid, or 
at aj = 0.9a (dashed curves), in which case it flows with a large, but finite viscosity. 



perhaps be directly tested experimentally by field-effect flow control of electro-osmosis in mi- 
crofluidic devices H58[|159lll60 ^, if a nearly constant double-layer voltage could be maintained 
along a channel and then varied systematically for diff'erent electrolytic solutions. Alternatively, 
one could compare with molecular dynamics simulations of electro-osmotic flow in nanochan- 
nels with oppositely and highly charged walls, leading to large induced double-layer voltages. 



A general MHS slip formula — If we assume s = s^'m Eq. ( 54 ) with a = = 1 (an arbitrary 
choice for analytical convenience), then for any MPB model with |p*^^| < oo we can integrate 



Eq. ( [49| ) obtain an MHS formula for the efl'ective electro-osmotic slip outside the double layer, 

^eff-'^o±—^-^p±^ (56) 

where ± = sign(^) = -sign(^z)) and pd^o) and qi^o) are given by ([17]) and ([18]), respectively. 
This simple and general expression reduces to the HS formula {^eff = ^d) in the limit of point- 
like ions (p^ax ~^ ^) and/or small voltages Q¥d ^ kT/e). At large voltages, \^eff \ is always 
reduced, depending on how the model for yu^^ imposes volume constraints or other contributions 
to the excess chemical potential. 

In any MPB model, such as Bikerman's, where counterions (±) form nearly uniform con- 



densed layers with c± ^ c^^^ and r] T]h from Eq. (53 ), the apparent zeta potential (56) either 
saturates to a constant, 

^.//~T* = -— In(^) (57) 

z+e \ Co I 
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if the condensed layer is stagnant (p'j = p^ax^ switches to a slower linear dependence 



a+ 



±\3 



(n-^c) 



(58) 



if the condensed layer has a finite viscosity, as shown in Fig. [17] The former case pi) ressembles 
the logarithmic concentration dependence of equilibrium zeta potential observed in many mi- 
crofluidic systems |262|, as well as the decay of ICEO flow noted above. The latter case allows 
for intermediate behavior between strong saturation of ^eff ~ the HS limit ^eff = ^d- 

MHS slip with Bikerman's model — The reduction of ^eff arises in diff'erent ways depending 
on the difl'use-layer model. In Bikerman's model for a symmetric electrolyte |57|, the limiting 
behavior is reached fairly suddenly. In that case the integral ( [55] ) can be performed analytically 
to obtain a simple formula: 



In 



1 + 4a Co sinh 



2kT 



(59) 



illustrated in Fig. [17] For a rigid condensed layer, this model predicts a simple logarithmic decay 
of ICEO flow with concentration, Eq. ( [57] ). If the condensed layer has a large, but finite viscosity, 
the decay is slower and more complex via Eq. (58 ). A general feature of even these very simple 
models is that ICEO flow becomes concentration-dependent and ion- specific at large voltages 
and/or high salt concentrations, through z±,a+, aj,and cq. 

It is interesting to compare Eq. ( [59] ) to the only previous MHS formula of Lyklema and Over- 
beek 163] [64l, based on the viscoelectric eff'ect ([52]) in the context of PB dilute-solution theory. 



As shown in Figures 14 and 17 the two formulae make similar predictions of saturation of the 



zeta potential with voltage with aj = a, but our formula ( 59 ) is simpler and more amenable to 
mathematical analysis, as illustrated below. (In contrast, the LO formula takes several lines to 
write down in closed form 1641 .) The parameter a is also more constrained on physical grounds, 
to be of order the hydrated ion size, than the empirical viscoelectric constant /. Qualitatively, 
the LO formula based on PB theory does not off'er any explanation of the experimental fact that 
ICEO flows depend on the particular ions, even in difl'erent solutions of the same ionic valences, 
{Zi} (such as NaCl and KCl), although we do not claim that our Bikerman MHS formula correctly 
captures any specific trends. 

MHS slip in a hard-sphere liquid — As noted above, hard-sphere liquid models show qual- 
itative diff'erences with Bikerman's lattice-based model, beyond just allowing the use of more 
realistic, smaller eff'ective ion sizes. As shown in Fig. [7] steric eff'ects are stronger in a hard- 
sphere liquid. The ion density approaches close packing as the voltage is increased much more 
slowly with the Carnahan- Starling model when compared to Bikerman's model, as shown in 
Id). When the CS model is used to compute the eff'ective zeta potential with pj = p*^^ and 
1 in Eq. ( [53] ) as shown in Fig. [18] the eff'ective zeta potential does not saturate as with 
Bikerman's model, and the layer continues to flow until extremely high voltages are reached, 
albeit much more slowly than in the classical HS model without charge-induced thickening. 

However, at this time we have no reason to assume that a = JS = 1 in Eq. ( [53] ). The crucial 
exponent controlling the viscosity divergence is J3. If we set = 4, then the divergence is fairly 
strong as shown by the dashed curves in Fig. [18] Here, we find a strong saturation of ^^ff^ similar 
to what is predicted by Bikerman's model (Fig. 17]), which more easily forms a condensed layer. 
A range of possible MHS slip behavior is possible with a given MPB model for the ion density 
profile, depending on the precise postulate for charge-induced thickening. 
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Figure 18: Effective zeta potential ^eff versus diffuse-layer voltage at different bulk concentrations using the 
Carnahan-Starling MPB model for charged hard spheres of diameter a = A Afrom Fig. [t] The concentration are 
Co = 1//M, 1 mM, and 1 M from top to bottom. The solid curves use the MHS slip formula {53) with a = /3 = I 
and pj = pmax, and the dashed curves change to j3 = 4. 



It is important to note that the CS equation was never intended to be used at high volume 
fractions and only fits the homogeneous chemical potential in the liquid state for O < 0.55. 
From Eq. [27]and Fig.jTJd) we see that the chemical potential diverges as O ^ 1, but the physical 



maximum for random close packing for bulk hard spheres is O 0.63 ll272M273ll271L which is 
exceeded in the CS MPB theory already at 1 volt across the diff'use layer. Moreover, crowding in 
the double layer occurs against a hard wall, which removes geometrical degrees of freedom and 
thus reduce the accessible local volume fraction for random close packing of hard spheres I293L 
and induces correlations not captured by the LDA | 250|. It is therefore not clear what the proper 
value of p*^^ should be to control our postulated viscosity divergence. The development of 
accurate models of the local rheology of highly charged double layers should ideally be guided 
by molecular theories and simulations. We simply give a range of examples to show what kind of 
qualitative behavior can be predicted by various simple MPB/MHS models, which are suitable 
for macroscopic theory and simulation of electrokinetic phenomena. 

Comparison to compact-layer models. For completeness, we briefly discuss how the general 
liquid-state models we develop above compare to the traditional approach of dividing the double 
layer into two parts, a flowing difl'use layer and a rigid compact layer, and assuming a constant 
viscosity and permittivity everywhere, leading to the HS formula. First we consider the classical 
Gouy-Chapman- Stern model which postulates an uncharged dielectric monolayer of solvent of 
constant efl'ective thickness As in contact with a difl'use layer of point-like ions obeying PB 
theory. The nonlinear charge- voltage relation I129L 

^ = ^, + 2^sinhf^V (60) 



zeAs \ 2kT 



then implies that only a logarithmically small diflfuse-layer voltage contributes to the zeta 
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Figure 19: Compact layer effects in traditional HS slip theory assuming a thin dielectric coating between the surface and 
the diffuse layer, whose importance is controlled by the parameter S = As I Ao{cq) = Cd{cq)ICs defined in Eq. |6}. Both 
PB (dot-dash curves) and Bikerman's MPB (solid curves) models are considered for the diffuse layer. The zeta potential, 
equal to the diffuse-layer voltage ^ = ^/), is plotted versus the total double layer voltage ^ at different values of the bulk 
salt concentration cq (labeled) for (5 = 0.1 in (a) and 6 = 10 in (b). 



potential at large total voltage ^ applied to the double layer, 

T 1 

(61) 
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This HS model also predicts that the zeta potential decays logarithmically with concentration and 
somewhat resembles our MHS models, as shown in Figure [19] 

Although the thin-dielectric compact-layer hypothesis leads to some reasonable predictions, it 
is not fully satisfactory either. As noted above, the saturation of zeta depends on the pileup of 
point-like ions, albeit reduced by transferring most of the voltage to the compact layer. This can 
be avoided by using an MPB model with steric constraints, such as Bikerman's model, together 



with a Stern layer and HS slip. As shown in Fig. 19 a), this tends to reduce the Stern-layer 
effect, since the diffuse layer is able to carry more voltage due to its reduced capacitance. With 
increasing S, however, the Stern layer always carries most of the voltage, and the difference 
between PB and MPB models on HS slip with the Stern model is eventually lost, as shown in 



Fig.[T9];b)for^= 10. 

This exercise shows that some concentration and voltage effects on slip can be captured by 
the classical Stem model and HS slip formula, but we are left with the same general criticisms 
made above for charging dynamics. The model places most of the large voltage across a region 
for which no detailed physical model is assumed, without attempting to account for liquid-state 
properties, salt concentration, surface roughness, liquid- surface interactions, etc. It is also not 
clear that a hypothetical Stern monolayer of solvent could withstand several Volts, while leaving 
a dilute diffuse layer with a small voltage, and some model for dynamical effects in response to 
applied voltages should be required. 

We have already seen in section [2] that the classical model is unable to predict all the trends in 
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nonlinear electrokinetics. Except in the case of a true dielectric coating (e.g a native oxide on a 
metal electrode), it seems more physically realistic to discard the classical concept of the com- 
pact layer (defined as the inner region where PB theory does not apply) and instead ascribe all 
of the double layer response to the dynamics of the liquid state, except perhaps for the dielectric 
response of a true Stern monolayer solvating the charged surface. An appropriate modified the- 
ory of the liquid can approximate the traditional properties of the compact layer through strong 



molecular interactions, as sketched in Fig. 15 At least in this work, we have shown that it is 
possible to describe a variety of nonlinear eff'ects in charging dynamics and electro-osmotic flow 
at large voltages and/or large salt concentrations without resorting to lumping the errors from 
dilute solution theory in a hypothetical compact layer outside the continuum model. Perhaps a 
more accurate theory would combine the ideas of this paper for the liquid phase with boundary 
conditions representing a compact interfacial phase, along the lines of the dynamical Stern layer 
model 12881 [2901 [29T1 [2921 . 

4.2.2. Ion-specific electrophoretic mobility 

We have already noted that induced-charge electrokinetic phenomena are sensitive to the so- 
lution composition in our models, via both the nonlinear capacitance and the efl'ective zeta po- 
tential. There are some surprising, general consequences, which are well illustrated by ICEP of 
an uncharged metallic sphere. It has been predicted using low- voltage models (H) [23l [24l and 
observed |294, 25 1 that asymmetric polarizable particles in a uniform field have an ICEP velocity 
scaling as E^, but it is widely believed that linear velocity scaling, U = bepEh, implies nonzero 
total charge. For non-polarizable particles with thin double layers, this is the case, unless the par- 
ticle has both asymmetric shape and a non-uniform charge density [|2951 12961 . For polarizable 
particles, our models predict that a nonzero mobility bgp can result simply from broken symmetry 
in the electrolyte, even for a perfectly symmetric particle. 

Consider an ideally polarizable, uncharged sphere of radius R in an asymmetric binary elec- 
trolyte with ions of unequal eff'ective sizes a+ and charges z+e, subject to a uniform DC back- 
ground field Eh. The first eff'ect to consider is the shift in the potential (po(Eh) of the particle 
(relative to the applied background potential) due to its asymmetric nonlinear capacitance, since 
the induced charge (which must integrate to zero) is more dense on one side than the other, as 



shown in Fig. 20 As noted above, this already yields a nonzero mobility with the HS formula, 

bep = Eb(f)o/Vb- 

The eff'ect of zeta saturation ( [57) ) provides a diff'erent dependence on the solution composition, 
which dominates in large fields, E^ » kT/eR. Since \(f)o(Eh)\ ^ E^R, the change in polarity 
of the induced charge occurs near the equator, around which there is only a narrow region with 
\^eff \ < l^c I- this limit, therefore, we approximate one hemisphere with uniform ^eff = 
and the other with which yields the ion-dependent mobility 
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The apparent zeta potential from a DC electrophoresis measurement ^ep = Vb^epl^, tends to a 
constant of order ^ kT/e, independent of E^ and R.lnsiZ'Z electrolyte, the limiting value 

3 kT ^ a+ 

In- (63) 

2 ze a- 



is set by the size ratio a+/a- and does not depend on the bulk concentration cq. With equal sizes 
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(a) (b) 

Figure 20: Mechanisms for DC electrophoretic motion U of an uncharged metallic sphere in an asymmetric electrolyte 
{62) due to saturated ICEO flow {57) in a large field Ejj » kT/eR. (a) Larger cations (below) in a z : z electrolyte 
(a+ > a-, z+ = \z-\) pack at lower density than smaller anions (above) and thus cover more of the sphere, but produce 
less slip due to greater crowding, (b) Divalent cations (below) cover less area and also produce less slip than monovalent 
anions of the same size (a+ = a-, z+ > \z-\). In both examples, the sphere has an apparent positive charge (U > 0). 



a+ = a- = a, the limiting apparent zeta potential 

IkT 

2 ze 



o — In(^^o) (64) 



is set by harmonic mean of the valences, z = z+Z-/(z+ + Z-), if z+ + Z- 0. 

An interesting feature of this nonlinearity is that the limiting mobility is set by properties of 



the electrolyte and is independent of the true charge of the particle. As sketched in Fig. |20| the 
induced viscosity increase alone causes the neutral sphere to have an apparent charge whose sign 
is that of the ions which condense at a lower potential (larger z and/or larger a). For consistency, 
however, we should also include the nonlinear capacitance effect, which can act in the opposite 
direction, making the apparent charge that of the ions which pack less densely (smaller z and/or 



smaller a). As shown in Fig.|2T] it turns out that the nonlinear capacitance effect is stronger and 
determines the sign of the apparent charge of the particle in large fields and/or high salt concen- 
trations. Nevertheless, the charge-induced thickening effect significantly reduces the mobility in 
this regime and introduces a strong decay with salt concentration. 

By now, it should be clear that ICEP of charged, asymmetric, polarizable particles can have 
a very complicated dependence on the solution chemistry in large electric fields and/or high salt 
concentrations. We must stress again that we do not include strongly nonlinear effects such as 
surface conduction and salt adsorption by the double layers, so the predictions of this paper only 
pertain to moderate voltages and thin double layers in the weakly nonlinear regime. Nevertheless, 
we already see some interesting new qualitative features predicted by the modified models. Our 
examples also show that the electrophoretic mobility of a homogeneous polarizable particle need 
not provide a reliable measure of its total charge, contrary to common wisdom. 
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Figure 21: Electrophoretic mobility of an ideally polarizable sphere of total charge Q in an asymmetric electrolyte with 
{a-la+f = 10 from Fig. [l3] with HS slip (dashed curves) compared to the same calculations redone with MHS slip jso) 
for (a) high salt concentration y_ = 0.1 and y+ = 0.01 and (b) moderate salt concentration y_ = 10~^ and y+ = 10^^^^. 
Charge-induced thickening has the opposite effect of nonlinear capacitance, since it gives more weight to smaller and/or 
more highly charged counterions in determining the electrophoretic mobility. This significantly reduces the apparent 
charge of the particle at large fields and/or high salt concentrations, but not enough to change its sign. 



4.2.3. Concentration dependence of AC electro-osmosis 

Next we revisit the weakly nonlinear analysis of ACEO pumping by adding the effect of 
charge-induced thickening via the MHS slip formula ([59]) in a symmetric electrolyte. In Fig. 



22 we show predictions of an ACEO pump including steric effects (Bikerman) and the simplest 



MHS, Equation 59 with ajla = 1. In interpreting these data, it is useful to remember that our 
simple model predicts the effective zeta potential to be the same as the voltage across the double 
layer up to a critical voltage whereafter the effective zeta potential saturates, see Fig. [T7] At high 
frequency, there is insufficient time for the double layers to fully charge and therefore the cutoff 
voltage is not reached and the slip is uneffected. At low frequency and high volatge, the double 
layers fully charge and the saturated zeta-potential severly limits the flow. Thus the MHS model 
acts as a high-pass filter for electroosmotic slip. This effect can be seen in Fig. 22 as the ion 
size decreases (thus the cutoff voltage increases). The upper dashed curve, given the physical 
interpretation of the model would require an ion size of 0. 1 angstroms at 0. 1 mM concentration. 
For more realistic ion sizes (at 0.1 mM and a = 3 angstroms, v = 1 x 10"^) we find that the 
viscoelectric effect essentially eliminates the prediction flow reversal from Bikerman's model. 
To date we have not successfully predicted flow decay with concentration and high frequency 
reversal in ACEO with a single, unified model though work continues in this direction. 

Figure [23] shows results applying the MHS model, ([59]), to the ACEO pumps of Urbanski et 
al f40 |. We see a decay in the maximum flow velocity with concentration that is reminiscent 
of experiments. Fig. [3] when we view the data in dimensionless form. Figure 23 (a) shows the 
dimensionless frequency response as we change concentration as in the experiments using a fixed 
ion size of 4 nm. As with the predictions of flow reversal in ACEO, it seems that the ion size that 
best fits the experiments is an order of magnitude larger than we would expect. The simulated 
data show the promise of a simple charge-induced thickening model to predict decay of flow with 
concentration. It is interesting to note that at the relatively low voltage of these experiments (3 
Vpp ^ 1 Vrms) there is no flow reversal predicted even when the viscoelectric effect is relatively 
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Figure 22: Frequency response of ACEO pumping with MHS slip and Bikerman's model of steric effects. Data are 
shown in dimensionless form where frequency is scaled by the RC time of the equivalent circuit and the velocity is 
scaled by sV^/t]Lo |42 |. The values of y = 2coa^ used in ^ are 0, IQ-^^, 10"^, IQ-"^, and lO'^, from top to bottom. 
Even at the extremely low values of y, this model of charge induced thickening essentially removes the prediction of 
reverse flow. For the steric efl'ects, v = 0.01, in the Bikerman model was used for all cases. 




Figure 23: Predicted frequency response of the ACEO pump of Urbanski et al f40] from Fig. |3ja) with an MPB 
double-layer model also accounting for MHS slip with charge-induced thickening. In (a) we show the dimensionless 
frequency response of the pump as we change concentration (C = 0.001, 0.003, 0.01, 0.03, 0.1 0.3, 1, and 10 mM from 
top to bottom). In (b) the same data are plotted in dimensional form. The model of charge-induced thickening uses a 
constant ion size of 4 nm. The response is computed with the Bikerman model using a constant value of y = 0.01, though 
no high-frequency flow reversal is predicted at this voltage. 
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weak, contrary to the experimental observations. 

We show the data in dimensional form in Figure 23 (b) for a direct comparison to experiments. 
We again see one of the key discrepancies which occurs in all models be they based in classical 
electrokinetic theory or based on the modified models presented here. All the models assuming 
blocking electrodes predict that the features in the frequency response are strongly concentration 
dependent, while the data show relative insensitivity to concentration. As discussed in Section 



3.2.1 we hypothesize this discrepancy is due to neglect of Faradaic reactions though this remains 



an area for future study. 

These simple models are capable of reproducing at least in qualitative way the general ex- 
perimental trends of ACEO. The Standard Model which does not account for charge induced 
thickening does not predict a strong flow decay as concentration is increased as all ACEO ex- 
periments have shown. We should also emphasize that our models are surely oversimplified, and 
various neglected eff'ects, such as specific solvent-mediated forces or electrostatic correlations 
can eff'ectively increase the range of crowding eff'ects (to allow the use of realistic ion sizes) and 
change their form in ways that may improve the ability to fit experimental data. 



5. Mathematical modeling of electrokinetics in a concentrated solution 

5.7. Nanoscale physics 

The fundamental difl&culty in modeling all electrokinetic phenomena is that complex 
molecular- scale phenomena at the electrified interface give rise to macroscopic fluid motion. 
The modified models above attempt to take into account some new physics - steric eff'ects for 
finite-sized ions and charge-dependent visco-electric eff'ects - but these are just simple first steps 
away from dilute-solution theory in electrokinetics. In this section, we develop a general theo- 
retical framework for modeling electrokinetic phenomena in concentrated solutions (including 
large applied voltages in dilute solutions). 

We begin by discussing various neglected eff'ects in our models above that might need to be 
incorporated into the general theory, some of which were anticipated (qualitatively) by J. J. Bik- 
erman decades ago 12971 12981 . In some sense, what we are attempting is to model part of the 
"compact layer" in microscopic detail using the same continuum equations that describe the 
"diff'use layer" and the bulk electrolyte. The hypothetical partitioning of the interface between 
diff'use and compact parts has become entrenched in electrokinetics and had many successes, 
at least in describing linear phenomena with non-polarizable surfaces of fixed charge. In such 
cases, it usually suffices to define an eff'ective slip plane, which marks the sharp transition from 
a dilute solution to a stagnant 12621 [TqI or dynamical 12881 12901 129 1 1 12921 compact layer, whose 
fixed position can be fit to electrokinetic measurements. At large induced voltages or concen- 
trated solutions, however, we believe it is necessary to describe the nanoscale rheology of the 
liquid in more detail, since it is otherwise not clear how to shift the eff'ective shear plane with 



voltage or concentration, as sketched in Fig. 15 As in prior work, it may still be useful to main- 



tain the theoretical construct of a separate "compact layer" via eff'ective boundary conditions on 
the continuum region, but crowding eff'ects, which vary with the local electric field and ionic 
concentrations, should also be included in modified electrokinetic equations. 

The following are some nanoscale physical eff'ects, other than volume constraints and vis- 
coelectric eff'ects (in the local density approximation), that we have neglected or included only 
heuristically, which may be important in nonlinear electrokinetics and other situations discussed 
below. In some cases, simple continuum models are available but have not yet been applied to 
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electrokinetic phenomena as part of a coherent theoretical framework. Setting the stage for such 
modehng is the goal of this section. 



• Electrostatic correlations. To our knowledge, all mathematical models in electrokinetics 
are based on the mean-field approximation, where each ion only feels an electrostatic force 
from the mean charge density of all the other ions In reality, ions are discrete charges that 
exert correlated forces on each other, which become especially important with increasing 
valence |187, 186|. The breakdown of the mean-field approximation for multivalent ions 
can lead to counterion condensation on the surface, eff'ectively leading to a correlated two- 
dimensional liquid (or glass) resembling a Wigner crystal of like charges in the limit of 



strong coupling ||279l [18911. 

In addition to this eff'ect, a promising direction for continuum modeling (discussed below) 
may be to build on MPB equation of Ref. I299L which describes the eff'ective restoring 
force acting on an ion that tries to fluctuate from its local electrostatic equilibrium position 
(in a one-component plasma 11891 ), only to be drawn back toward its correlation hole. 
We conjecture that this efl'ect may not be so important at very high charge densities in large 
applied voltages, where simple crowding becomes dominant and flattens out any oscillations 
in the charge-density profiles. Some evidence comes from recent molecular simulations of 
the double-layer capacitance of ion liquids 18311 , which verify the square-root scaling of 
MPB theory with steric eff'ects discuss above. 

On the other hand, electrostatic correlations may be crucial for the dynamics of highly 
charged double layers under mechanical shear stress, which to our knowledge has never 
been studied. We have conjectured above that the correlation hole interaction may eff'ec- 
tively enhance the viscosity of the solution and lead to charge-induced thickening. Electro- 
static correlations may also eff'ectively increase the critical length for crowding eff'ects, and 
we have noted elsewhere that the Bjerrum length £b = e^/4nsiykT is at the same scale as 
the eff'ective ion size required in our steric MPB models of ACEO, especially if corrected 
for reduced permittivity in large field 1601 11301 [67l . The electrostatic correlation length Ac 
(defined below) is approximately Ac = z^^b for a z : z electrolyte, so ion-ion correlations are 
particularly important for multivalent ions. 

The relative importance of corrections to mean-field theory due to electrostatic ion-ion cor- 
relations is controlled by the dimensionless parameter, 

= 17 = Ms,kT)y^ electrolyte) (65) 

which grows with bulk salt concentration like ^/co, as the Debye screening length shrinks. 
The Bjerrum length in bulk water is £b ~ 7 A, so we expect strong correlation eff'ects 
on electrochemical transport and electrokinetics when A^ > 7z^ Awhich is 7 A, 1.4 nm, 
and 2.8 nm for monovalent, divalent, and trivalent ions, respectively. The condition of 
"intermediate coupling" Sc = 0(1) is met in many concentrated aqueous solutions, so it 
may be necessary to include correlation eff'ects in theories of nonlinear electrokinetics and 
double-layer charging dynamics. Below, we discuss a possible modification of Poisson's 
equation |299 J, which could provide a starting point. 
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• Specific ion-ion interactions. In addition to entropic effects such as hard-sphere repul- 
sion and long-range electrostatic forces, ions can interact via more complicated short-range 
forces, due to direct molecular interactions or solvent-mediated effective forces. DiCaprio 
et al. 1 193] recently expressed the excess free energy density of an electrolyte as the sum 
of a hard-sphere entropic contribution (e.g. the CS model above, or its generalization to 
polydisperse mixtures) plus the first, quadratic terms in a Taylor expansion in the ionic 
concentrations. In our theoretical framework, the latter corresponds to an additional linear, 
enthalpic contribution to the excess chemical potential 

11^^ = ^ UijCj (specific ion-ion interactions) (66) 

j 

where the coupling coefficients Uij (with units of volume) are assumed to be concentration 
independent. These terms effectively "soften" the hard- sphere interactions by introducing 
further short-range forces. For example, if an > (or < 0), then like-charged ions ex- 
perience an additional repulsion (or attraction) within a volume an, and this increases (or 
decreases) the effective hard- sphere radius for crowded counterions. (More complicated 
concentration-dependent enthalpic terms have also been postulated for ions intercalated in 
crystalline solids in rechargeable battery electrodes |300|. ) Specific interactions may also 
contribute to dynamical friction coefficients between different ionic species, e.g. leading to 
off-diagonal elements in the mobility tensor (see below), which could be important in large 
ac fields 1 176] where oppositely charged ions must quickly pass each other upon polariza- 
tion reversal. 

• Ion-surface correlations. Just as ions in the liquid phase can have short-range interactions 
with each other, beyond the mean-field electrostatic force and hard-sphere repulsion, so too 
can they have specific interactions with molecules comprising a phase boundary, such as 
a solid wall. The simplest of these result from hard- sphere ordering and solvent-mediated 
forces which produce correlations near the flat wall. Molecular simulations taking into 
account the solvent (beyond the primitive model) reveal strong density oscillations near a 
surface 13011 11051 11061 , and we have already discussed how non-local WDA theories can 
capture hard-sphere aspects 12441 1245| 1195) [2491 1196^ Solvent-induced lay- 
ering can also be described by local continuum MPB models, by simply adding an effective 
external potential V^i^) to the excess chemical potential I106II1941|108L 

^f^ = V„{x) = -kT\niP^)^ (67) 

which can be expressed in terms of the statistical density profile (pair correlation function) 
psix) of the solvent molecules near the wall, relative to their bulk density [3011. Contri- 
butions to the effective ion- wall potential V^ix) can also arise from other effects, such as the 
polarizability of ions 13021 and "cavitation forces" due to disruption of hydrogen bonding 
networks in water 13031 1304 1: these effects have been studied for water/gas interfaces, but 
could also be important at solid metal surfaces. 

Stern's original model of a solvation monolayer separating ions from the surface 11831 can 
be viewed the simplest model of type ( [67] ) to account for excluded volume, since he essen- 
tially postulated an infinite chemical potential V^vix) for ions in the monolayer. Of course, 
it is more realistic to allow for smooth oscillations in the ion- wall correlation function over 
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several molecular diameters, but at least Stern's model is easily incorporated into a boundary 
condition 1 129. 1541 (see below). 



• Surface heterogeneity. Discrete surface charges contribute to electrostatic correla- 
tions I279L but chemical heterogeneities can also play an important role in the structure 
of the double layer and electrokinetic phenomena I3Q51 13061 . In the context of PB theory, 
the effect of surface roughness on the differential capacitance of the double layer depends 
on the correlation length of the roughess relative to the Debye length B071 13081 . With 
finite-size ions, roughness could have a much stronger effect, not only on capacitance, but 
also on slip generation. Molecular dynamics simulations of electro-osmosis over atomically 
rough surfaces have revealed departures from PB theory 12801 . but large induced voltages 
have yet to be studied. Given the arguments of the previous section, the charge-induced 
viscosity increase is likely to be enhanced by roughness, since the electrostatic compression 
of hydrated ions against molecular scale asperities should thicken the fluid and make the 
apparent ion size seem larger. 

• Specific adsorption of ions. Another effect we have neglected is the specific adsorption of 
ions, which break free of solvation and come into direct molecular contact with a surface, 
as already included by Stern 11831 in his model for the isolated equilibrium diffuse layer 
based on a Langmuir adsorption isotherm. Much more refined surface adsorption models 
have been developed since then, especially to describe the colloid chemistry of oxidic ma- 
terials 13091 . Because the specific adsorbed of ions, and thus the surface charge, is a direct 
function of the ion concentration within the diffuse layer directly next to the interface, the 
close approach of two colloidal particles and the overlap of their respective diffuse layers 
will lead to surface charge modulation, so-called "charge regulation" 13101 13051 . For iden- 
tical interfaces approaching one another, the surface charge will be reduced, and in the limit 
of touching Stern (adsorption) planes, the surface charge will become zero. Charge regula- 
tion of different interfaces (hetero-interaction) leads to the sum of the two interfacial charge 
densities to approach zero when the interaction distance is reduced, an effect which can be 
so strong that for amphoteric materials (i.e., those that can be both positive and negatively 
charged) it will lead to the inversion of the surface charge on one of the interfaces. Simulta- 
neously, the force-distance curve can be highly non-monotonic with for instance repulsion 
at contact and at sufficient separation, but with an electrostatic attraction at intermediate 
separation BIOL 

Ion adsorption effects are begining to be considered in electrokinetics and may be partic- 
ularly important in nonlinear electrokinetics, due to the large driving force for desolva- 
tion in large applied voltages. Effects of specific ion adsorption on the equilibrium surface 
charge of metal electrodes as function of surface potential Bill and on the electro-osmotic 
flow in a microchannel between electrodes have been considered by Duval |312|. Adsorp- 
tion is an important mechanism of ion specificity of electrokinetics with hydrophobic sur- 
faces B13II304|| . In the present context of induced-charge electrokinetics, we have already 
highlighted the recent work of Suh and Kang 11751 11791 incorporating surface adsorption 
of ions in models of ACEO flow. 

In electrochemistry, the effect of specific ion adsorption is widely invoked to explain the 
increase of differential capacitance of the double layer at high voltage observed in many 
experimental situations, since the distance between plates of the equivalent capacitor effec- 



tively shrinks from the "outer Helmholtz plane" to the "inner Helmholtz plane" (Fig. 24), 
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Figure 24: Sketch of the double layer near a blocking electrode at high voltage. Solvated counterions (green) are 
crowded in the inner region and smoothly transition across the outer diffuse region to a dilute solution with solvated 
anions (orange). An ion can break free from its solvation shell and adsorb on the surface (black), thereby moving from 
the outer Helmholtz plane (OHP) to the inner Helmholtz plane (IHP) 



and a number phenomenological models are available H55[|314ll . Some results combining 
specific adsorption with MPB models of the diff'use layer are in Ref. I26QI . Of course, at 
very large voltages with steric eff'ects included, the diff'erential capacitance at blocking elec- 
trodes must eventually decrease, once the IHP and OHP both become saturated with ions, 
as we have argued above. This regime of universal square-root decay is often inaccessible 
at low frequency due to Faradaic reactions, but our modeling of ACEO flow reversal above 
suggests that it can be probed at high frequency. 

• Normal current and Faradaic reactions. Faradaic electron-transfer reactions can afl'ect 
local surface potentials and ion concentrations, and thus also electro-osmotic flows over 
electrodes. This efl'ect has recently been studied by Duval et al. I312II3151I3T6,.3171 in the 
context of nonlinear streaming potentials over electrodes in the presence of a redox couple 
in solution. Reversible redox couples, in particular, can greatly reduce streaming potentials 
over electrodes due to conduction that results from bipolar electrolysis |316|, and this efl'ect 
can be isolated using an indifl'erent supporting electrolyte II317L These studies are based on 
the standard Butler- Volmer equation, with reactions driven by the full double layer voltage 
(rather than just the compact part | 318 |) and a Frumkin correction based on dilute-solution 
theory, so some modifications may be required for large applied voltages (see below). 

Most models for induced-charge electrokinetics assume blocking surfaces in order to focus 
on capacitive double-layer charging and simple ICEO flow. As noted in section 2, however, 
there is growing evidence that Faradaic reactions play a major role, especially at low fre- 
quency and high voltage (as in Fig. [TT]). This leads to normal currents, which can perturb the 
equilibrium structure of the double layer. Lacoste et al. 1 22 1 have recently noted that reverse 
ICEO flows can arise even at low voltages in a mean-field theory of biological membranes 
passing a normal ionic current, in the Helmholtz limit of a thick dielectric layer (6 oo). 
Since this limit corresponds to a large "correction factor" in the Standard Model as inferred 
from most experiments (see section 2), it may be that non-equilibrium double-layer struc- 
ture in the presence of Faradaic reactions is involved in the low-frequency, high- voltage 
flow reversal in ACEO and TWEO, in addition to efl'ects of difl'usion-layer electroconvec- 
tion L151I . Below, we note that the mathematical description of Faradaic reactions may also 
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need to be modified for large voltages or concentrated solutions. 

In the next section, we summarize a general mathematical framework for the dynamics of elec- 
trolytes and ionic liquids, which naturally allows the incorporation of some of these eff'ects. 

5.2. Modified electrokinetic equations 
5.2.1. Continuum modeling approaches 

Until now we have focused on situations with thin double layers and integrated over the 
double-layer structure to obtain modified eff'ective boundary conditions on the quasi-neutral bulk 
fluid, namely the MPB difl'erential capacitance in section [3] and MHS slip in section |4] In this 
section, we summarize the modified continuum electrokinetic equations corresponding to these 
models, which could be applied to arbitrary geometries with thin or thick double layers in non- 
linear electrokinetics. A few examples of this approach been developed for linear electrokinet- 
ics L68*i69J and electrochemical dynamics 1130112211 . Such modified continuum equations aim 
to capture more of the essential physics of nanoscale electrokinetics, while remaining much sim- 
pler and widely appHcable than brute-force molecular dynamics ifTM 11051 [TQ6|[T071[TQ81[TQ91 
or more complex statistical approaches 1319113201177] . and could have broad applicability beyond 
the problems considered here. 

We have made the case for modified electrokinetic equations based on experimental and the- 
oretical arguments in nonlinear electrokinetics, but similar issues also arise in other fields. We 
have already noted that Bikerman's model has recently been adapted to model the double layer 
in ionic liquids |[8T| [82l |83l, but mainly to predict the diff'erential capacitance for use in RC 
circuit models. We are not aware of any attempt to describe electrokinetics or non-equiilibrium 
dynamics of ionic liquids, so our general approach below may have some relevance, in the limit 
of very high salt concentration, approaching a molten salt. Steric eff'ects in polyelectrolytes have 
also been described by Bikerman's model |220| as well as the CS and BMCSL hard-sphere 
models 12351 12341 , and this is another interesting area to consider electrokinetics with crowding 
eff'ects. 

Even in the more familiar context of electrolytes, there has been a recent explosion of interest 
in nanofluidics |[89l|93|, since the high surface-to-volume ratio of nanochannels amplifies the 
importance of transport phenomena occurring in confined geometries, eff'ectively inside double 
layer. In recent years, the classical electrokinetic equations for a dilute solution have been used 
extensively (and exclusively) to model nanochannel phenomena such as ion selectivity 191] [92| 
[93l[94l and mechanical-to-electrical power conversion |!95l l96l[97l[98l[99lfTQQl , but quantitative 
agreement with experiments often requires fitting "extra" compact-layer properties. Modified 
electrokinetic equations with additional physics might improve theoretical predictions without 
relying as much on adjusting boundary conditions. 

As noted in the introduction, Liu et el. [77 1 have recently considered correlations and crowd- 
ing eff'ects in steady nanochannel transport via a continuum MPB theory, based on more com- 
plex statistical thermodynamical models than those considered here ^TSj |79l [801 • This approach 
incorporates more physics than we have considered above, but this comes at the expense of sacri- 
ficing some generality and mathematical simplicity, since it assumes equilibrium charge profiles 
and requires numerical integration of coupled integral equations to determine the self-consistent 
charge density in Poisson's equation. Nonlinear integral equations for ion profiles or correlation 
functions also result from other statistical approaches II241L such as the hyper-netted chain and 
related models |186, 187], and statistical density functional theory 1 3191 12391 13211 . For hard- 
sphere models, perhaps the simplest theories of this type are based on the weighted-density ap- 
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proximation f250l, as discussed above. In spite of many successes, however, all of these methods 
require significant eff'ort to solve numerically, even in simple situations, while their underlying 
physics remains relatively simple (e.g. charged hard spheres). As such, it may be more fruitful 
to apply modern computers to molecular dynamics simulations with more realistic interatomic 
forces than to solve non-local continuum models numerically. 

Twenty years ago, in the context of fitting double-layer capacitance data, Macdonald con- 
cluded that "integral-equation statistics treatments are too complicated and too limited... to be of 
practical usefulness", while lattice-gas models may be "entirely adequate to describe the diff'use 
layer" 11881 . Although computational advances have made integral-equation approaches more 
feasible today, we believe that simple mathematical descriptions are still useful, if not necessary, 
to model non-equilibrium dynamical phenomena. Therefore, we focus here on local continuum 
models with finite-size ions, following Cervera et al. 1681 [69l . We will also mention a simple 
non-local WDA, which could serve as a first correction to the LDA in the general theory. 



5.2.2. Electrochemical transport 

We now present a general modeling framework based on electrochemical potentials, which 
applies to non-equilibrium situations and includes the simple cases considered above. In princi- 
ple, one could start with the full theory of non-equilibrium thermodynamics of multicomponent 
systems 1322. 323 L but we develop a simpler phenomenological theory for electrokinetics in an 
isothermal concentrated solution. Our general starting point is to postulate a MF-LDA contin- 
uum model for the electrochemical potential of an ion of species / (possibly including a solvation 
shell, depending on the model), decomposed into ideal, electrostatic, and excess contributions as 
follows: 

lii = kT In Ci -\- ZiC^ -\- iif{x, {Cj}, 0) (68) 

The thermodynamic meaning of jdi is the Gibbs free energy diff'erence upon adding a particle 
of species / (and replacing other particles or empty space) within a continuum element, viewed 
as a local open system in quasi-equilibrium with the reservoir of nearby elements. Its gradient 
-VyU/ acts as a mean "thermodynamic force" on each particle, driving the system toward local 



equilibrium. In (68), we have defined jdf as the excess chemical potential of ion / relative to 
a dilute solution, expressed in terms of local continuum variables, such as the position x (e.g. 
distance from a surface), ionic concentrations [cj] and their gradients (to approximate non-local 
contributions), electrostatic potential and field E = -V0, etc. The fundamental significance 
of electrochemical potentials is emphasized by Newman |190], who also questions validity of 
the mean electrostatic potential cp at the molecular level in a multi-component concentrated so- 
lution. Nevertheless, it is necessary to separate long-range electrostatic forces from short-range 
chemical interactions to develop continuum equations for electrokinetics, so we proceed with the 



phenomenological decomposition in (68). From this perspective, we can view yu^^ = kTXnfi as 
defining the chemical activity coefficient / in terms of the mean-field approximation for 0. 

In principle, the excess chemical potential jd^^ can be derived from microscopic statistical mod- 
els or by fitting to molecular dynamics simulations or experiments. For example, throughout this 
article, we have focused on two simple MF-LDA models of ji^^ for steric effects of excluded vol- 



ume for solvated ions, namely Bikerman's lattice-based model (20) and the Carnahan- Starling 
hard- sphere-liquid model ( [27] ). Above in this section, we noted some other possible contribu- 
tions. Simple MF-LDA models are also available for specific solvent-mediated ion- wall interac- 
tions ( [67] ) and ion-ion interactions ( [66] ). 
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More complicated effects of statistical correlations can also be included in //^^ by making it 
a non-local functional of the ion density-profiles. The simplest prescription of this type is the 
weighted density approximation (WDA), where jd^^ for a bulk homogeneous liquid is evaluated 
with local reference concentrations in place of c/, which are obtained by non-local averaging 
of the nearby inhomogeneous concentrations, 



Several methods to construct the weight function w(r) are available for general liquids F 2441 12451 
12461 12481 1 189 J and have been applied in equilibrium MF-WDA electrolyte simulations 112491 
11951 11961 11971 125QI . For hard- sphere liquids, the simplest and most natural choice for w{r) 
is the hard- sphere indicator function, w = 3/(47ra^) for r < Ui and w = for r > a/, which 
turns out to be quite successful, so it would be interesting to try to include this particular WDA 
in electrokinetic models, as a first correction to the LDA. This WDA leads to realistic density 
oscillations near a surface without the need to fit an empirical wall potential ([67]) in the LDA. On 
the other hand, LDA models are much simpler to implement numerically and allow analytical 
progress, while any non-local continuum model may be intractable for dynamical problems or 
complicated geometries. 

In non-equilibrium thermodynamics 13221 13231 1 190L the mass flux densities F/ are obtained 
from the phenomenological hypothesis of linear response: 



where = -CjVjdj is the thermodynamic force density (force/volume) acting on species j and Lij 
is the (symmetric, positive definite) Onsager mobility tensor converting these forces into mean 
drift velocities in the frame moving with the mass-averaged velocity u. The mobility tensor is 
related to the diff'usivity tensor by the Einstein relation, Dij = LijkT, and is usually assumed 
to be be diagonal, Dij = DiSij, but this can only be justified for a dilute solution. In a highly 
concentrated solution, the mobility tensor (or its inverse, the friction tensor) may have significant 
off'-diagonal elements I3221l323ll . 

The Onsager mobility in ( [70| ) is generally take to be a constant (i.e. linear response), but in 
highly concentrated and dissipative liquids, there may be nonlinear concentration dependence. 
Various models for glassy relaxation 1275112741 and granular flow 12761 have postulated a power- 
law decay of the mobility. 



as the particle volume fraction O approaches the jamming or glass transition O^, respectively. 
The crowding and compression of counterions against a highly charged surface by a large, time- 
dependent normal electric field may cause a similar, temporary decrease in mobility very close to 
the surface. Such a nonlinear eff'ect on ion transport would be consistent with the charge-induced 
viscosity increase proposed above. 

With the fluxes defined by ( [70| ), the diff'erential form of mass conservation is 




(69) 




(70) 




(71) 



(72) 



where is the reaction rate density for production (or removal) of ion /, which is usually (but not 
always |257 |) set to zero for electrolytes. With these further assumptions, the modified Nemst- 
Planck equations take the general form II13QL 



Zie(p + jdl 
kT 



(73) 



5.2.3. Electrostatics 

The system of modified PNP equations I13QI is usually closed by making the mean-field ap- 
proximation, in which the electrostatic potential self-consistently satisfies Poisson's equation 

-V '(sV(P)=p = YjZieci (74) 

where p is the mean charge density. The permittivity £ of a polar solvent like water is usually 
taken as a constant in ([74]), but numerous models exist for field-dependent permittivity 6:(|V0|), 
such as ([34]) or ( [35] ), discussed in section |3.1.6[ The classical eff'ect of dielectric saturation 
reduces the permittivity at large fields due to the alignment of solvent dipoles 11551 12561 12281 
I188L although an increase in dipole density near a surface may have the opposite eff'ect 12541 . 
The permittivity can also vary with temperature, due to Joule heating or reactions, but here we 
only consider isothermal systems. 



We are not aware of any attempts to go beyond the mean-field approximation ( [74] ) in dynam- 
ical problems of ion transport or electrokinetics. This would seem to require a simple contin- 
uum treatment of correlation eff'ects, ideally leading to a general modification of ( [74] ). Recently, 
Santangelo 12991 derived a simple modified PB equation accounting for ion-ion electrostatic 
correlations in a one-component plasma 1 189 | near a charged wall in the relevant regime of "in- 
termediate coupling", Sc = Ac/ Ad = 0(1), which suggests modifying Poisson's equation with an 
additional term, 

{aIv^-1)V-(sV(P)=p (75) 

where Ac is the electrostatic ion-ion correlation length, set by the balance of thermal energy and 
Coulomb energy in the dielectric medium. Physically, the extra term roughly accounts for inter- 
actions between an ion and its correlation hole during thermal fluctuations. The higher derivative 
of the correction term introduces the possibility of oscillations in the ion densities at the scale of 
the correlation length. (It also requires additional boundary conditions, discussed below.) The 



relative importance of the correction term in ( 75 ) is measured by the dimensionless parameter 



Sc = Ac /Ad introduced above, which takes the form Sc = z^b/Ad for a z : z electrolyte. Since 
Sc = 0(1) for concentrated aqueous solutions (and increases if local permittivity decreases), cor- 
relation efl'ects could be important in electrokinetics, and Santangelo 's equation may provide a 
useful starting point for analysis. For ionic liquids, correlation efl'ects are even more important, 
since the difl'use layer shrinks to the molecular scale (Sc » 1). 

5.2.4. Electrochemical hydrodynamics 

To determine the mass-averaged solution velocity, u, we enforce the conservation of linear 
momentum 

Pm^ + V.T = f (76) 
ot 

where is the total solution mass density (which can be set to its bulk value in most cases), 
T is the hydrodynamic stress tensor (arising from mechanical friction), and f = 2; ft is the 
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thermodynamic force density (acting on all the ions, independent of fluid flow). The nonlinear 
inertial convection term u • Vu could be added to ( [76| ), but it is typically negligible in nonlinear 
electrokinetics due to a very small Reynolds number. For the stress tensor, the first approximation 
is the Newtonian form, T = pi - T^^^ with contributions from the dynamic pressure p (to satisfy 
incompressibility, V • u = 0) and the viscous stress tensor. 



\ dxj dxi 



(11) 



As explained in section |4j there may be significant local changes in viscosity near a highly 
charged surface. For example, the following empirical form combines the viscoelectric eff'ect 



(52) and charge-induced thickening ([53 



J / 



2V 



(78) 



For simplicity, one could typically set 7 = 0, since E tends to grow with p in a similar way, via 
the MPB equations, as noted above. Alternatively, one could set = 0, but that removes any 
explicit dependence on the eflfective ion sizes or other modified interactions. 



The thermodynamic force f, which acts as a source of momentum in (76), can be simplified, 
if we assume small departures from local thermal equilibrium. We also neglect heat transfer and 
assume isothermal conditions. In that case, the Gibbs free energy density, g = Yji Ciiih varies as 



(79) 



where po is the hydrostatic pressure. Taking the variation between adjacent continuum elements, 
we obtain 

-f=^c,V//,-V/7o+pV0 (80) 



which is a form of the Gibbs-Duhem relation 13221 11901 , adapted for an isothermal charged 
system. Note that po includes the osmotic pressure that balances concentration gradients, which 
takes the form, kTYuiCi, in a dilute solution upon inserting ([T3| into (80). In a concentrated 
solution, the osmotic pressure can take a more complicated, possibly non-algebraic form, but its 
gradient should still uphold the local Gibbs-Duhem relati on([80| near thermal equilibrium. 

Since we assume incompressible flow, we can insert (|80|) into (76) and absorb /7o into the 
dynamical pressure p. In this way, we arrive at the familiar form of the unsteady Stokes equation. 



(81) 



with an electrostatic force density, = -pV0. The unsteady term dxxldt in (81 ) is often over- 
looked, but it can be important in nonlinear electrokinetics, e.g. for oscillating momentum bound- 
ary layers and vortex shedding in response to AC forcing |[T6ll43]| . 

The Stokes equation ( [8T] ) is the standard description of fluid mechanics at low Reynolds num- 
ber, which is normally applied in to a dilute solution, but we see that it also holds for an isother- 
mal, concentrated solution near equilibrium, regardless of the form of yu^^. In their theory of 



non-equilbrium thermodynamics, De Groot and Mazur [3221 instead assert (81) as the funda- 
mental expression of momentum conservation, where the "external" or "long-ranged" force 
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acts as the source of momentum flux in a continuum element, whose internal "short-ranged" 



forces are described by T. Here, we show the equivalence of starting with (76) based on the full 



thermodynamic force f for a concentrated solution and deriving (81 ) as the quasi-equilibrium 
limit, where only the external (electrostatic) force produces momentum. In this limit, which is 
consistent with the assumption of linear response ( [7Q| ) for the mass fluxes, all other the "chemi- 
cal" interactions in f - only contribute to the (osmotic) pressure. 

Non-equilibrium thermodynamics can be extended to account for the electrical polarizability 
of a concentrated solution |322|. For a linear dielectric medium with variable permittivity the 
electrostatic force density can be expressed in the familiar form, = -V • T^^\ where 

T^^ ^s[EiE^-^-\E\H^ (82) 

is the Maxwell stress tensor [324]. As noted above, in polar solvents, the permittivity should gen- 
erally decrease in large fields. Various phenomenological models for e(E) can be incorporated 
into the theory of electrokinetics for concentrated solutions, but they complicate analysis and can 
introduce seemingly unphysical oscillations or singularities in the concentration profiles 12601 
and are perhaps best avoided, or included only heuristically in the boundary conditions. 

53. Modified boundary conditions 
5.3.1. Electrostatic boundary conditions 

For Poisson's equation in the mean-field approximation ( [74] ), the electrostatic boundary con- 
ditions at a dielectric surface require continuity of the tangential electric field and equate the 
jump in normal electric displacement across the interface to the free charge qs (which is re- 
lated to the equilibrium zeta potential) II324II . For a low-dielectric surface with a fixed surface 
charge density, the internal electric field can often be neglected, yielding the standard boundary 
condition, 

sh-V(f) = qs. (83) 

Alternatively, for a metal surface, one can simply fix the potential = 0o or allow for a thin 
dielectric layer (or compact Stem layer) on the surface through the mixed boundary condi- 
tion [iT54l 1771 fT44l [221. 

A(Ps =(P-(Po = ^sn-^(p-^, (84) 

where As = shs /ss is an eff'ective thickness of the layer, equal to the true thickness hs by the 
ratio of permittivities of the solution s and the layer ss , and Cs = £s/hs is its capacitance. The 
boundary condition can also be generalized for voltage-dependent surface capacitance, which 
makes the surface-layer voltage drop A^^ a nonlinear function of the normal electric field 11541 . 
When applying ([84]) to a metal electrode, one can set qs = to model the Stern layer as a thin 



dielectric coating of solvent molecules I228L while specific adsorption of ions would lead to 
qs ^0. 



The preceding boundary conditions can be imposed on the mean-field Poisson equation (74), 
but the modified equation for electrostatic correlations ( [75] ) introduces a fourth derivative term 
and requires one more boundary condition on each surface. Charge conservation requires the 
following boundary condition 

h-[{AlV^-l)sV(p]=qs (85) 

where brackets indicate the jump across the boundary. For consistency with the derivation of 
( [75] ), Santangelo sets the second term (jump in mean dielectric displacement) to zero and thus 
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equates the surface charge to the jump in the curvature of the field 12991 . For an insulating 
surface of fixed charge density, this would imply replacing ( 83 ) with two boundary conditions 



h-V^ = Q and h- A^V^ (sVcp) = qs (86) 
For a metal surface with a compact Stern layer modeled as a thin dielectric coating (with a 



uniform electric field), we would replace ( 84 ) with 



(p-(po = Ash-V(p and h • A^V^ (sV(p) = qs . (87) 

This mathematical model provides an interesting opportunity for analysis of correlation eff'ects 
in electrochemical dynamics and electrokinetics, although it is only a first approximation. 

5.3.2. Electrochemical boundary conditions 

Standard boundary conditions for the concentration fields equate normal ionic fluxes with 
surface reaction rates 

/T-F,=/?,({Q},0,{yU,},...), (88) 

which vanish for inert ions {Ri = 0). The flux boundary condition can also be generalized for 
a dynamical Stern layer, which supports tangential ionic fluxes 12881 1290112911 or adsorption of 
ions LI8OIII 81 1117511179 L although some efl'ects of this type are already captured by the modified 
electrokinetic equations, as noted above. The reaction rate Ri may describe surface adsorption, 
in which case there is an auxiliary equation for the surface concentration, 

^+V,.F,=7?,-, (89) 
ot 

where the second term allows for surface diffiision. If the kinetics of this reaction are fast (large 
Damkoller number) and surface transport is slow (small Dukhin-Bikerman number) compared to 
bulk transport, then the reaction is in quasi-equilibrium. In that case, the surface concentration 
Cs is given by an adsorption isotherm, which equates the surface and nearby liquid chemical 
potentials, [is = j^i. For example, the popular Langmuir isotherm follows from a lattice-gas model 
of the surface adsorption sites, yu^ = kT ln[Cs/(Cmax - Cs)], with dilute solution theory for yU/. The 
reaction rate Rf may also describe Faradaic electron-transfer reactions, such as electrodeposition 
(which also involves an adsorption step, and resulting motion of the metal surface) or redox 
reactions (which alter the charge of ions remaining in the liquid region). 

The proper mathematical description of electrochemical reaction kinetics is complex and not 
fully understood |155l . The standard model in electrochemistry is the Butler- Volmer equation, 
usually applied across the entire double layer under conditions of electroneutrality 11901 . Ap- 
plying an analogous expression at the molecular level has better theoretical justification and in- 
troduces the "Frumkin correction" for diff'use-layer voltage variations 13251 . See Ref. | 318 | for 
a recent review. For example, for the redox reaction R <r^ O -\- ne~, this model asserts Arrhenius 
kinetics for the forward (anodic) and backward (cathodic) reaction rates, 

R = kaCRC-''^'''^^''^^ - kccoe''^'''^'^''^^, (90) 

where ka and kc are rate constants for the anodic and cathodic reactions, cr and cq are concen- 
trations of species R and O, and aR and ao are transfer coefficients defined below {aR -h a^o = 1) 
. The bias voltage A^^ can be interpreted as the Stern-layer voltage in models of the type we 
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have considered here f l91[ 13261 13271 13181 . This approach has been used to model ACEO P2ll 
and TWEO 1 151 1 at reacting electrode arrays in dilute solutions. It is straightforward to include 
nonlinear differential capacitance of the Stem layer, Cs{/^(ps)^ as well 1191 L but more significant 
modifications may be needed for concentrated solutions and large voltages. 

For consistency with our theoretical framework based on non-equilibrium thermodynamics, 



the reaction rate should properly be expressed in terms of the electrochemical potentials 13281 . 

R = kQ [e^^-^^Ts)lkT _ ^(Mo-M'^ykT^ , (91) 

where jir and jdo are the complete electrochemical potentials of the reaction complex in the 
reduced and oxidized states, and yu^^ is the excess electrochemical potential in the transition 
state, and is an arbitrary rate constant (which can be set by shifting yu^^). The Butler- Volmer 



equation (90) follows from dilute- solution theory (//^^ = /j,^ = 0) and a purely electrostatic 
model for the activation barrier which is a linear combination of the electrostatic energy of the 
reduced and oxidized states, weighted by the transfer coefficients: 

= Ea-\- aRqRCp -h ao(qo(p - necpo) (92) 

where Ea is a composition-independent activation energy barrier, absorbed into the rate constants 
ka and kc, and qR and qo = ^r-^ ne are the charges of the reduced and oxided states. The general 



expression ( [91] ) can be derived from statistical transition- state theory in a hypothetical local open 
system, and it contains a variety of possible non-electrostatic influences on the reaction rate, via 
the excess contributions to the chemical potentials. This approach was recently introduced in 
the context of ion intercalation in rechargeable-battery materials, where the chemical potential 
in the electrode, and thus the reaction rate, depends on gradients in the ion concentration [ 300] . 
Steric efl'ects were also recently considered in the context of fuel-cell membranes I327L and a 
form of ( [9T] ) was efl'ectively applied. In nonlinear electrokinetics, it may also be necessary to 



consider more general forms of the reaction rate in (91 ), whenever the voltage is large enough to 



invalidate the dilute- solution approximation close to the surface. 



5.3.3. Hydrodynamic boundary conditions 

Until recently, almost all theoretical studies in electrokinetics have assumed the no-slip bound- 
ary condition for the liquid velocity, u = U, where U is the velocity of the surface. With the emer- 
gence of microfluidics [TJ, the phenomenon of hydrodynamic slip has been studied extensively 
in simple, Newtownian fluids 12631 12641 12661 and interpreted in terms of the Navier boundary 
condition (3291, 

Au = u-U = /?/t- Vu, (93) 

where the slip Au is proportional to the shear strain rate via the slip length b. Flow past smooth 
hydrophilic surfaces has been shown to be consistent with the no-slip hypothesis, but b can reach 
tens of nanometres for hydrophobic surfaces [267', '268', '2691 or even several microns over super- 
hydrophobic textured surfaces with trapped nanobubbles 13 301 1331 113321 1333] [334 1. 

The study of electrokinetic phenomena in the presence of slip was perhaps first pursued by 
the group of N. V. Churaev 13351 [336 1. For electro-osmotic flow in a microchannel, Kiseleva 
et al. [335.1 considered the eff'ect of exponentially varying viscosity ri{x) near a wall, increasing 
toward a hydrophilic surface or decreasing toward a hydrophobic surface, and Muller et al. [13361 
studied the impact of the slip boundary condition ( [93] ), which enhances the flow by a factor 
(1 -h blAi)) at low voltage [12641 . This enhancement of electro-osmosis was recently analyzed and 
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demonstrated by Joly et al. fl081 via molecular dynamics simulations and extended to dilfusio- 
osmosis by Ajdari and Bocquet |337|. (The permittivity may also vary near a hydrophobic 
surface, and this also can affect particle interactions II338II .) The possibility of slip-enhanced 
(linear) electro-osmotic flows has generated considerable excitement in nanofluidics 11001 13391 , 
but so far it has only been analyzed with the classical electrokinetic equations and the simple, 
purely viscous boundary condition ( [93] ). 

We suggest using a modified Navier slip boundary condition I340L 



Au = M(T-/2), (94) 

where T - his the total normal traction on the surface due to short-range forces (force/area) and 
M is an interfacial mobility tensor (velocity x area/force), which is non-diagonal for anisotropic 
surfaces. For an isotropic, impermeable surface, the mobility matrix is diagonal, M = MI, with 



zero elements for normal flow, and then Eq. (93 ) is recovered from (94 ) with b = Mrj. 



Using these models, it would be interesting to study the competition of viscoelectric efl'ects 



[78]) and hydrodynamic slip ( 93 ) or ( 94 ) in nonlinear electrokinetic phenomena at polarizable. 



hydrophobic surfaces in large applied voltages and/or concentrated solutions. The simple low- 



voltage amplification factor (1 -h b/Aoico)) associated with (93 ) increases with concentration and 
becomes appreciable when the diff'use-layer thickness becomes smaller than the slip length, so 
hydrophobic surfaces may counteract the eff'ect of increasing viscosity very close to the surface 
and allow induced-charge electro-kinetic phenomena to be observed at higher concentrations and 
voltages. On the other hand, these predictions depend sensitively on the location of the slip plane 
where the hydrodynamic boundary condition is imposed, and how it relates to the compact-layer 
plane (OHP, Stern plane, reaction plane, etc.) where the electrochemical boundary conditions 
are imposed, especially in the concentrated- solutions models. For example, we have seen that 
charge-induced thickening reduces the flow in the difl'use layer, but this can be counteracted by 
hydrodynamic slip if the slip plane lies closer to the surface than the thickened region, due to the 
amplification of viscous stress on the slip plane. We pose the eff'ect of hydrodynamic slip in a 
highly charged double layer on a metal surface as an interesting open question for future work. 

5.4. Thin double layers and diffusion layers 

The modified electrokinetic equations and boundary conditions above may be useful in model- 
ing nanoscale electrokinetic phenomena, e.g. taking into account steric eff'ects of finite ion sizes, 
but at larger scales, where the double layers become thin, matched asymptotic expansions can 
be used to systematically integrate out the diff'use layer and derive eff'ective boundary conditions 
on the quasineutral bulk. First, we briefly summarize the results in the typical situation where 
the voltage is not strong enough to drive the double layer out of equilibrium or fully deplete the 
bulk salt concentration, due to difl'usion limitation. In that case, the ion transport equations ( [72| ) 



remain unchanged in the bulk, but Pois son's equation ( [74| ) is replaced by the condition of elec- 
troneutrality, Y^iZieCi = 0. The fluid equations are also unchanged, and bulk viscosity variations 
can usually be neglected. 

As noted above, chemical potentials are approximately constant (or "quasi-equilibrium" holds) 
in the normal direction across a thin double layer. For the ionic concentrations, the boundary 
conditions then take the form of surface conservation laws II114L 

^+V,-F«=«-F,-/?,, (95) 
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where Ti(^D, {ci}) is the excess concentration of species / per unit area, V^-F[^^(^z), {c/}) is the sur- 
face divergence of the integrated tangential flux in the difl'use layer, h • F/ is the normal flux from 
the bulk, and Rii^o, {ci}) is the Faradaic reaction rate density at the surface, evaluated in terms of 
the bulk variables, which thus includes the Frumkin correction. See Ref. 1 1 SI for expressions for 
Ti and F|^^ using Bikerman's model, neglecting convective fluxes, and Refs. 1 1541 13261 134 11 13 181 
for expressions for Ri for Faradaic reactions in a dilute solution, neglecting tangential transport. 
Equation ( [95] ) generalizes the RC boundary condition ([2]). 



Integrating over the difl'use-layer also yields a "first kind" eff'ective slip boundary condition 
for the bulk fluid velocity. 



u, = b^'""^ E, - y b^!^'^ kTVt In ^ , (96) 
where the first term describes electro-osmosis driven by the bulk tangential field with 



b^^^\^^D,{ci}) given by (49) and various approximations above. The second term describes 



diff'usio-osmosis in response to tangential bulk salt concentration gradients [ 3421 11121 11771 [H 
fT23]| . The diffusio-osmotic mobilities bf''\^D, {ci}) can be systematically derived from the full 
transport equations above in the limit of thin double layers, following the asymptotic analysis of 
Prieve et al. II177II or Rubinstein and Zaltzman I.123>,124J . Simple expressions can be derived for 
dilute solution theory, such as 



= — 
rib 



<E,-4||fl„c„sh(|£)v,l„i. 



(97) 



for a dilute z ' z electrolyte, where c = c+ = c_ is the quasi-neutral bulk salt concentration. 
(Equivalent forms are in Refs. I1771ll23ll . and this corrects a missing factor of 4 in Eq. (Ill) of 
Ref. |115 |.). More cumbersome expressions, which are not expressible in terms of elementary 
functions, result from modified equations with volume constraints I26QI . Similar asymptotic 
methods can also express the eff'ective slip boundary condition in terms of bulk electrochemical 
potential gradients, {V^yU/}. 

In our many examples of induced-charge electrokinetic phenomena with blocking surfaces, 
we have assumed "weakly nonlinear" dynamics 11291 1131L where the bulk concentration is 
not significantly perturbed. Under "strongly nonlinear" dynamics at large voltages, even with 
blocking electrodes 111291 13431 13441 11741 11761 11751 1179L strong bulk concentration gradients 



can develop, and the other terms in Equations ([95]) and ([96]) become important I I 12l 11231 112411 . 



Although steric eff'ects generally reduce the importance of surface conduction (smaller Dukhin- 
Bikerman number) compared to dilute solution theory since there are not nearly as many ions in 
the double layer (601, diff'usio-osmosis, concentration polarization, and diff'usion-layer electro- 
convection 1 1511 are aff'ected less and could be significant. In addition to Faradaic reactions, 
other sources of normal ion flux, such as salt adsorption by the difl'use-charge layers and Lang- 
muir adsorption of ions on the surface, can also produce time-depedent difl'usion layers (thicker 
than the difl'use charge layer, but thinner than the bulk region) oscillating at twice the frequency 
of the AC forcing, although these layers are weakly charged and may not have a large efl'ect on 
the efl'ective fluid slip lUSl 11221 \lM- All of the problems we have considered above should 
be revisited in the strongly nonlinear regime to better understand the predictions of the modified 
electrokinetic equations, but this is beyond the scope of the paper. 

As noted in section 2, the eff'ects of Faradaic reactions or other mechanisms for normal ionic 
flux are still poorly understand in nonlinear electrokinetics, even for thin double layers. Normal 
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currents can disturb the quasi-equilibrium structure of the double layer, even at small currents 
(in the Helmholtz limit S 0), and lead to seemingly reverse ICEO flows |22|. Concentration 
gradients can also develop due to normal ionic fluxes under diff'usion limitation. If the bulk salt 
concentration approaches zero, at a limiting current, the quasi-equilibrium difl'use layer expands 
into a non-equilibrium space charge layer and drives second-kind electro-osmotic flows |113J 
and hydrodynamic instability II1221 1 12411 . The classical electrokinetic equations suflftce in that 
case, since the decreasing concentration only helps to validate dilute solution theory. In the 
case of induced-charge electrokinetic phenomena, however, the strongly nonlinear regime is just 
beginning to be explored and may require modeling with modified electrokinetic equations. 



6. Conclusion 

We have provided a critical review of recent work in nonlinear "induced-charge" electroki- 
netics, comparing theory to experiment for the first time across a wide range of phenomena to 
extract general trends. In doing so, we were naturally led to question the theoretical foundations 
of the field and develop modified equations for electrokinetics in concentrated solutions at large 
voltages. These equations may find applications in diverse areas of electrochemistry and fluid 
mechanics. A crucial aspect of this efl'ort was to survey microscopic models of electrolytes from 
difl'erent fields, where convection and charge relaxation are neglected. Our experience shows the 
importance of integrating knowledge across scientific communities. For example, we managed 
to find seven independent formulations of Bikerman's model (1942) from 1947 to 1997 and two 
rediscoveries of Freise's capacitance formula (1952) in 2007, each leading to separate literatures 



without any cross references (section 3.1.2). 



We have argued that many new nonlinear phenomena can arise at large induced voltages, and 
we focused on three that could play a major role in induced-charge electrokinetics: (i) Crowding 
of counterions against a blocking surface (Fig.|4]) decreases the diff'erential capacitance (Fig. [5]), 



which may explain high frequency flow reversal in ACEO pumps (Fig. 10) and imply ion 



specific mobility of polarizable particles in large fields (Fig. [12]); (ii) a charge-induced viscosity 



increase upon ion crowding (Fig. [15]) reduces the eff'ective zeta potential (Fig. [17]), which implies 
flow decay with increasing concentration and an additional source of ion- specificity (Fig. 20); 
(iii) Each of these eff'ects is enhanced by dielectric saturation of the solution in large electric fields 
(Fig. [9]). To illustrate these phenomena, we have derived extensive analytical formulae based on 
simple models in the mean-field and local-density approximations (MF-LDA), including lattice- 
gas and hard-sphere models for steric eff'ects of finite ion sizes, as well as various postulates of 
charge-induced thickening leading to modified electro-osmotic slip formulae. 

We have also developed a theoretical framework for electrokinetics based on non-equilibrium 
thermodynamics in a concentrated solution (section [5]). Although motivated by induced-charge 
electrokinetics, these general equations and boundary conditions could find applications in many 
other areas involving nanoscale electrochemical transport. Our examples of MF-LDA models 
focusing on the eff'ects above should be refined and extended in future work, e.g. to account for 
various solvent eff'ects, specific adsorption of ions, and Faradaic reactions, and we have reviewed 
some of the relevant literature for guidance. Especially for multivalent ions, it may be necessary 
to go beyond the MF approximation to account for electrostatic correlations, which we alter 
the inner double-layer structure and could contribute to charge-induced thickening. A proper 
description of volume constraints may require going beyond the LDA, especially very close to a 
surface, to account for short-range correlations and related density oscillations. 
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Legend has it that Wolfgang PauH once said,"if God made materials, then surfaces are the 
work of the Devil". At end of this study, it is tempting to draw the same conclusion, given 
the complexity of possible phenomena occurring at large surface potentials and our inability 
to devise a single model to predict all the experimental data. We remain optimistic, however, 
that simple, predictive models will follow from improved nanoscale understanding of the double 
layer. 

Our study raises an important general question for continuum modeling, "Where should the 
continuum region end and give way to a boundary condition?". This question has received much 
attention recently in the context of hydrodynamic slip, but the situation is much more compli- 
cated for electrochemical relaxation and electrokinetic phenomena. The pioneering papers of 
Stern 1 183 J and Bikerman |59 | introduced two opposing, general perspectives, which provide 
the historical context for our work. Stem first proposed describing the outer "diffuse part" of 
the double layer with dilute solution theory, while lumping any discrepancies into an empirical 
boundary condition on the inner "compact layer". Bikerman then showed that a similar finite-size 
cutoff* of dilute-solution theory could instead by obtained by modifying the bulk equations for a 
concentrated solution, without ajdusting the boundary condition. By the latter half of the century. 
Stern's approach became widely adopted through simple empirical models for the compact layer, 
and Bikerman's paper was essentially forgotten. 

We have argued that Bikerman's perspective should be revisited, especially for extreme con- 
ditions of large applied voltage, salt concentration, and/or frequency, since otherwise it is not 
clear how compact-layer boundary conditions should change to describe a dynamical region of 
nonlinear response near the surface. Bikerman's perspective also applies more easily to non- 
trivial nanoscale geometries, where the classical concepts of "Stern plane" and "slip plane" are 
less well defined. We have shown that simple continuum equations for ion crowding, dielectric 
saturation, and viscoelectric response all predict that the compact layer and slip plane both ef- 
fectively advance into the solution with increasing surface charge, but without the need to define 
these empirical concepts. On the other hand. Stem's perspective remains attractive for mathemat- 
ical modeling, since it preserves simple equations for the ffuid domain and lumps complicated 
molecular details into boundary conditions, so perhaps a combined approach is needed. Cer- 
tainly, true surface eff'ects, such as specific adsorption and Faradaic reactions, require eff'ective 
boundary conditions, consistent with any modifications of the bulk equations (section |5]2]). 

In our opinion, it remains a grand challenge to describe double-layer structure, electrochemi- 
cal relaxation, and ICEO ffow at a polarizable surface over vast parameter ranges, from 10 mV 
to 10 V in voltage, to 100 kHz in frequency, and 1 jjM to 1 M in ionic strength, using a simple 
- but not over-simplified - mathematical model, amenable to analytical results (in idealized lim- 
its) and numerical simulations (for typical experimental situations). The upper extremes of these 
conditions correspond to a new regime for the theory of electrokinetic phenomena, where counte- 
rions become crowded during time-dependent relaxation and ffow near a highly charged surface. 
Nanoscale experiments and molecular dynamics simulations will be crucial to further develop 
and validate various theoretical postulates. If properly validated, we believe that modiffed math- 
ematical models can advance our understanding and better predict experimental observations. 
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